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Abstract 



In this work a systematic method for determining correlated wavefunctions of ex- 
tended systems in the ground state as well as in excited states is presented. It allows 
to fully exploit the power of quantum-chemical programs designed for correlation 
calculations of finite molecules. 

Using localized Hartree-Fock (HF) orbitals (both occupied and virtual ones), an 
effective Hamiltonian which can easily be transferred from finite to infinite systems is 
built up. Correlation corrections to the matrix elements of the effective Hamiltonian 
are derived from clusters using an incremental scheme. To treat the correlation 
effects, multireference configuration interaction (MRCI) calculations with singly and 
doubly excited configurations (SD) are performed. This way one is able to generate 
both valence and conduction bands where all correlation effects in the excited states 
as well as in the ground state of the system are taken into account. An appropriate 
size-extensivity correction to the MRCI(SD) correlation energies is developed which 
takes into account the open-shell character of the excited states. 

This approach is applicable to a wide range of polymers and crystals. In the 
present work trans-polyacetylene is chosen as a test system. The corresponding 
band structure is obtained with the correlation of all electrons in the system being 
included on a very high level of sophistication. The account of correlation effects 
leads to substantial shifts of the " center-of-mass" positions of the bands (valence 
bands are shifted upwards and conduction bands downwards) and a flattening of all 
bands compared to the corresponding HF band structure. The method reaches the 
quantum-chemical level of accuracy. 

Further an extention of the above approach to excitons (optical excitations) in 
crystals is developed which allows to use standard quantum-chemical methods to 
describe the electron-hole pairs and to finally obtain excitonic bands. 
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Chapter 1 
Introduction 



The study of electronic properties of solids was constantly a hot topic in physics 
over many decades. Various different approximations, e.g. the free-electrons model, 
the tight-binding method and others (see, e.g. |Z72j ). were tried to get qualitative 
description of the electron states in periodic systems. More sophisticated methods 
aiming already at a quantitative description have been developed since computers 
have become a common tool of researchers. Among these methods one could dis- 
tinguish a number of approaches based on the local density approximation (LDA) 
to density functional theory (DFT) |HK64j . |KS65j which are widely used presently 
in solid-state physics. They allow, in particular, to account for electron correlation 
effects by locally treating electrons as a homogeneous gas. These methods have 
become very popular as they provide quantitatively correct results for the ground 
state of many systems. However, these methods fail to describe correctly the band 
structure of solids, e.g. they usually give too narrow band gaps for non-conducting 
systems. To overcome this problem, further improvements of the formalism have 
been developed. For instance in the last decade the most successive methods became 
those based on the so-called GW approximation |H65j . 

One of main advantages of methods based on DFT is their relatively small com- 
putational cost that opens the possibility to routinely apply them to periodic systems 
with large unit cells. However, all DFT-based methods are focused on computing 
physical ground-state properties of the systems which are related to electron density 
and avoid calculation of the many-body wavefunction which would provide us with 
explicite information on the electron correlation in particular materials. Another 
serious disadvantage of DFT methods is the impossibility to control the accuracy 
of the obtained results and to improve them in a systematic manner. Thus, in the 
long perspective, an alternative way to study electron states in periodic systems is 
required. 

In parallel to the development of DFT-based approaches, methods which solve 
the Hartree-Fock (HF) equation for solids were established. They provide single- 
particle wavefunctions and energies obtained self-consistently with the many-body 
wavefunction being a superposition of Bloch waves. For this purpose powerful pro- 
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gram packages like CRYSTAL |PDR88j , |CRY98j have been designed. The use of the 



variational method to solve the HF equation allows one to reach the desired accuracy 
within a given set of parameters (such as basis sets, crystal geometry, thresholds for 
integral values etc.) and a systematic improvement of the basis sets and a reducing 
of the thresholds leads to further improvement of the results. However, in this ap- 
proach electrons are considered to experience only a mean interaction field and thus 
the electron correlation is missing totally. As a result of this, the fundamental band 
gap of non-conducting materials is substantially overestimated by the HF method. 

Let us now turn to quantum chemistry. In this field, an intensive development of 
methods which start with the HF solution for finite systems and accurately account 
for electron correlation took place. These methods yield correlated wavefunctions 
and correlated energies of finite molecules in good agreement with experimental 
data. A detailed description of such approaches can be found in, e.g. |S096j and 
|LM86j . Among these methods one can distinguish the configuration interaction 
(CI) methods, the coupled cluster (CC) methods and methods based on the many- 
body perturbation theory (MBPT) such as the second-order M0ller-Plesset (MP2) 
method which have become standard in quantum chemistry and are implemented 
in quantum-chemical program packages like MOLPRO |M2000j and GAUSSIAN 



At first glance, there is no possibility to apply these quantum-chemical methods 
(which operate with single-electron wavefunctions defined in real space) to extended 
systems where the states are usually defined in momentum space. However, if one 
reformulates the HF solution in terms of localized Wannier functions, which are 
defined in the real space, and considers finite clusters of the crystal, then one can 
account for correlation effect within this cluster. An electron placed in a localized 
orbital interacts strongly with electrons in neighbor localized orbitals trying to avoid 
coming too close to each other. Electrons in farther orbitals are already screened 
and "felt" via the mean field established in the HF description. Thus, electron 
correlation is predominantly a local effect, and therefore the electron correlation 
effects in an entire crystal can be derived from finite clusters. Here, very accurate 
quantum-chemical methods can be employed. One only needs a scheme to transfer 
the results from clusters to periodic systems which in fact resembles a Wannier 
back transformation. This way one gets the correlated wavefunction of extended 
system with which any physical quantity of a system can be calculated. The scheme 
described above preserves the spirit of wavef unction-has ed methods which consist 
of a sequence of well- controlled approximations and thus allowing to estimate and 
improve the accuracy of the final results. 

The implementation of this idea has started in the late eighties of the last century 
and from that time has reached a valuable success. Ground-state and cohesive 
energies of many covalent and ionic systems have been calculated until now (see e.g. 
Ht92i|, [Et92EI, |PF^95j . jPF^96j . jDDF^97j . jKPF^97j . jA^PnOaj and jA^DnObj ). 



The approach was also used to study structural properties of polymers f |YKD97] and 
IAA+991 ). In parallel to the study of ground-state properties of extended systems. 
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an approach for the hole states in periodic infinite systems has been developed. To 
the present moment it was addressed to the valence bands of simple semiconductors, 
such as diamond, silicon and germanium f |GSF93j . |GSF97j and |AFS00j ). Further 
information on the current state of the art of wavefunction-based approaches for 
extended systems which employ quantum-chemical methods can be found in the 
recent review |F02j . 

Another promising wavefunction-based method appeared recently and has suc- 
cessively been applied to the band structure of LiF (see |AI01j . |AF02j and |A02j ). 
In the post-HF step it uses Green's function formalism which is usually not consid- 
ered as a standard quantum-chemical method. Therefore, we only mention it here 
briefiy. 

As one can notice, the wavefunction-based methods for solids described in |GSF93j . 
|GSF97j and |AFS00j study correlation of electrons in occupied orbitals. However, 
to also obtain correlated conduction bands one needs a proper description of cor- 
relation effects for states with an additional electron in the virtual space as well. 
Without this one can not evaluate one of the most fundamental electronic property 
of insulators and semiconductors, namely the band gap. The main reason of the 
previous approaches to fail here is the impossibility to construct localized virtual or- 
bitals in clusters (in contrast to occupied orbitals) by standard localization schemes 
used in quantum chemistry. Therefore, a different approach to localization of virtual 
orbitals was needed. 

Conceptual works |MV97j and |SMV01j on obtaining maximally localized or- 
bitals in periodic systems and implementations of similar ideas in CRYSTAL code 
f lBZ+Ol] and |ZWDSOT| ) convinced us to use localized virtual orbitals obtained in 
a periodic system as a starting point for the construction of localized virtual molec- 
ular orbitals for sufficiently large cluster models of periodic systems. This way we 
could properly describe states with an attached electron in crystals by means of 
finite clusters. Thus, we can calculate not only valence but also conduction bands 
and, therefore, the band gap provided electron correlation effects are accurately and 
systematically taken into account. This we regard as the main motivation for the 
present work. 

In fact, having obtained the correlated hand structure we extract from it not only 
the band gap but also band widths, the ionization potential and the electron affinity 
of the system with electron correlation being taken into account. The correlated 
wavefunction is used implicitely to calculate these quantities. Our approach also 
allows us to estimate the error bar for each calculated quantity and it provides a 
way to improve the accuracy further. 

All the quantities calculated here can also be measured directly in experiments. 
To study the band structure of crystals photoionization (Fig. II. 1|) and inverse pho- 
toionization experiments fFig. II. 2p can be performed (see e.g. |PS99j ). 

Ionization potentials IP^ (the valence-band energies) are the difference between 
the energy of (A^ — l)-electron system with a hole in band n and the energy of the 
neutral A'^-electron system. This quantity can be measured by irradiation of light 
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Figure 1.1: A schematic view of photoionization process. 

with known frequency on the surface of a crystal and measuring the kinetic energy 
of the photoemitted electrons: E'kin = hu + {Eq — '^~^En) = hu — IP„. 

The inverse process is used to measure the electron affinities EA^ which corre- 
spond to the energy which a neutral system gains when an extra electron is accepted 
in the conduction band m and an (A^ + l)-electron state is formed. Electrons with 
known kinetic energy are shined onto the surface and the frequency of the emitted 
light is measured: = -Ekin + (-Eo ~ ^"'"^-E'm) = -E-kin + EAm- This two experiments 
give band energies which we aim to calculate theoretically. 

Another motivation of this study is to extend the borders of the wavefunction- 
based approaches to correlated band structures. By this we mean that we would like 
i) to investigate a system with more complex unit cell than was done before, ii) to use 
larger basis sets (valence triple zeta) which are necessary in quantum chemistry to 
obtain indeed accurate results and which are appropriate for correlation calculation 
in (A^ + l)-electron systems and iii) to use one of the most accurate and powerful 
correlation method (the multireference configuration interaction method restricted 
to singly and doubly excited configurations) available in standard quantum-chemical 
program packages. The implementation of all the ideas mentioned above put our 
approach far ahead with respect to other methods for band structure calculations 
in non-conducting systems. 

Yet another goal of our study is to develop an approach which allows to even 
describe well-bound excitons in crystals. These electron-hole pairs made up by neu- 
tral excitations of the electronic structure are rather compact objects and therefore 
can be treated in real space. Occupied and virtual Wannier orbitals are used for 
that purpose such that the electron-hole pair can be described in finite clusters of 
the crystal by standard quantum-chemical methods. Provided, the localized orbitals 
in such a cluster represent the Wannier orbitals of the crystal well, an appropriate 
transformation allows to transfer the relevant data on the localized excitonic states 
from the clusters to the periodic system and excitonic bands can be calculated. 

As a test system we have chosen trans-polyacetylene (tPA). This is a semicon- 
ducting polymer which consists of weakly interacting single chains. Therefore, its 
highly anisotropic electronic properties can be studied in a one-dimensional system, 
a tPA single chain. This is a very useful feature for the implementation of our 
approach to conduction bands in particular because the corresponding localized vir- 



Figure 1.2: A schematic view of inverse photoionization process. 



tual one-particle orbitals are in fact rather diffuse objects. To study the decay of 
the correlation effects with the distance of the involved localized orbitals we need a 
linearly scaling system. 

Also polymers themselves are intriguing systems providing researchers with sur- 
prising results. The discovery and development of electrically conductive polymers 
were appreciated by the Nobel Prize in Chemistry for 2000. One of the experi- 
mentally well-studied systems is doped trans-polyacetylene. Thus, the sound un- 
derstanding of the electronic structure of tPA is highly desired. In our work we 
also aim to answer the question whether electron correlation in polymers can be 
described properly in a single chain to give results which agree with experimental 
data for the bulk system. For this purpose, in section 3.1.1, we describe in details 
the investigated material and further, in section 4.4, we compare results obtained 
for tPA single chains with existing experimental data on bulk tPA crystals. 

The thesis is organized as follows. In Chapter 2 the basic theoretical aspects of 
our approach are discussed. At the beginning of the chapter we briefly sketch the 
approximations used in ab initio methods. In first sections we consider the solution 
of the HF equation for infinite systems in terms of Bloch waves and Wannier orbitals 
and define the local quantities which also show up in infinite systems and clusters and 
which can be transferred from one to another. In the last two sections we discuss why 
electron correlations are necessary to be taken into account, the physical meaning of 
these effects in both neutral and charged systems, and how these correlation effects 
can be evaluated by the MRCI(SD) method. 

Chapter 3 is dedicated to the development of our wavefunction-based approach. 
First we introduce the structural and electronic properties of the test system. Next 
we discuss in details how to get localized HF occupied and virtual orbitals in clusters 
and define the corresponding local matrix elements (LMEs) which are used to cal- 
culate the SCF band structure of tPA. In section 3.3 we discuss how to evaluate the 
correlation corrections to the LMEs by MRCI(SD) method, and the use of an ap- 
proximate but well controlled scheme for this purpose (namely incremental scheme) 
is presented in section 3.4. The possibility to reduce the size-extensivity error of 
CI(SD) is explored in section 3.5. In the last section of Chapter 3 the way to handle 
the problem of the emergence of satellite states in cluster models is presented. 

In Chapter 4 we summarize and discuss numerical results for ^rans-polyacetylene. 



6 



CHAPTER 1. INTRODUCTION 



There we emphasize the effect of electron correlation on the key quantities of our 
method, namely the local matrix elements. We use correlation corrected LMEs to 
plot the correlated band structure from which we get the band gap, IPs, EAs, and 
band widths of linear tPA. The accuracy of the whole method is explored and the 
error bars for the calculated quantities are estimated. In section 4.4 we compare 
our results obtained for the one-dimensional system to experimental data on bulk 
irans-polyacetylene. 

In Chapter 5 we introduce a new approach to excitonic states in crystals in terms 
of singly excited configurations. The method is established and some preliminary 
results are presented. 

In Chapter 6 we conclude and discuss perspectives of our approach. 



Chapter 2 

Local approach to wavefunctions 



In this chapter we would hke to discuss the basic theoretical approaches to our study. 
The ultimate goal of the physicist studying theoretically electronic properties of an 
atom, a molecule or even a bulk material would be to get a many-body wavefunc- 
tion which maximally approaches the exact one and contains all the information 
about the electronic structure, the forces on the nuclei, and other parameters of 
the system in all its possible states. Unfortunately, obtaining the exact many-body 
wavefunction is unrealistic. In order to proceed one has either to make some ap- 
proximations for the wavefunctions and then to calculate the desired quantities with 
this approximated wavefunctions (and to estimate the emerging errors) or one tries 
to obtain those quantities directly (avoiding the calculation of wavefunctions) using 
some tricks or semiempirical methods. In the latter case the question about the 
achievable accuracy often remains open. When one is interested not only in a quali- 
tative description of a system and the effects in it, but also in accurate quantitative 
results then the first way seems to be more attractive, even though more compli- 
cated. Having chosen this way, one choses the most accurate among all feasible 
approximations in order to get a realistic description of the system and to reduce 
the errors to a minimum. 

In our study we are aiming at a corre/atec? wavefunction of a periodic system (e.g. 
trans-polyacetylene). This means that we introduce some controlled approximations 
at the very beginning, then get the wavefunction which is assumed to be very close to 
the real one, and finally calculate the ground state energy and energies of one-particle 
states (the band structure). By the term "controlled approximation" we mean that 
we can both estimate the error bar for the calculated quantities and improve our 
results to (in principle) any desired accuracy by increasing the computational efforts. 
Here we mention these approximations. 

The first one is the Born-Oppenheimer approximation. The system under con- 
sideration consists of nuclei vibrating around their equilibrium positions and elec- 
trons filling the space between the nuclei. The Hamiltonian of the system and, 
correspondingly, its wavefunction depends on coordinates of both all nuclei and 
all electrons. As electrons are much lighter than nuclei they move much faster and. 



7 



8 



CHAPTER 2. LOCAL APPROACH TO WAVEFUNCTIONS 



therefore, can be considered as moving in the field of fixed nuclei (for details we reffer 
to |Z72j . chapter 6.11). Then solving only the "electronic" part of the Hamiltonian 
one obtains electron eigenfunctions and eigenvalues which depend parametrically on 
the nuclear coordinates. Having the solution for the electron problem one considers 
the nuclei as slowly moving in the averaged field of the fast moving electrons and 
solves the problem of vibration, rotation and translation of the system consistently 
with the help of electronic solutions. In our work we are concentrating on the elec- 
tronic problem and consider the nuclei as fixed at their equilibrium positions which 
we can define e.g. by using experimental data on the structure of a given mate- 
rial. By this we exclude the vibrational and rotational degrees of freedom of finite 
systems baring in mind, however, the possibility to obtain them by looking at the 
electronic structure for different positions of the nuclei. This approximation is valid 
as long as electron-phonon interaction in the system is negligible. 

The second point to discuss is the use of basis sets for the electron wavefunctions 
which are necessary for numerical calculations of electronic structures. We use 
atomic basis sets since quantum chemistry has demonstrated for many decades that 
for atoms and molecules this representation gives very accurate results as compared 
to experimental data. We also have to notice that this basis sets are finite. To set 
up an analytic expression for the wavefunction of an electron we use Gaussian-type 
atomic orbitals (GTO) which are very favourable from the computational point of 
view. With these GTO basis sets very accurate results can be obtained and the 
accuracy depends on the number of atomic orbitals in the basis set for each atom 
and on the number of Gaussians per atomic orbital. Enlarging these numbers leads 
to tremendous computational efforts when one wants to treat big molecules which in 
fact limits the size of basis sets. However, at this point we want to emphasize that 
finite but sufficiently large GTO basis sets can provide reasonable accuracy and are 
still convenient to use. 

In our work we consider systems consisting of light atoms whose atomic numbers 
are not larger than six. Therefore, we can neglect relativistic effects and consider the 
spin coordinates of electrons as independent from their spatial ones. Spin-orbital 
coupling only becomes not negligible for more heavy atoms starting from the atomic 
number around 30 (see, e.g. §72 in |LL77j and chapter 3.11 in |Z72j ). In zero-oder 



approximation we consider the wavefunction of a system as being the one of an 
independent particles system where electrons occupy spin orbitals (a function of 
electron position vector and spin coordinate at) or alternatively we can say that 
spatial orbitals (which depend only on rj) are occupied by pairs of electrons one 
having spin up and the other having spin down. Such spin orbitals are solutions 
of a single-particle Hamiltonian in the self-consisted field (SCF) approximation and 
are represented as linear combinations of atomic orbitals (LCAO). Then, the many- 
body wavefunction of the whole system is given as the antisymmetrized product of 
all occupied spin orbitals, the so-called Slater determinant (see e.g. (Ennsi). The 
ansatz to write molecular orbitals as LCAO is also an approximation and is based 
on the fact that when an electron is placed close to one nucleus its interaction energy 
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with this nucleus is much larger than interaction with all other nuclei. In quantum 
chemistry the LCAO approximation is very common and considered as being proved 
empirically to yield very good agreement between experimental data and theoretical 
results for molecules, once electron correlations are included. 

Despite the essentially unavoidable approximations mentioned so far we do not 
make any approximation for the nuclear potentials and the Coulomb interaction 
between electrons. The systems are treated as realistic as possible keeping by this 
the sense of ab initio calculations. 

As mentioned above we start from the wavefunction of the closed-shell system 
obtained in the SCF (or HF) approximation. In this approximation the one-particle 
Hamiltonian is constructed, treating an electron as moving in the field of the fixed 
nuclei and the averaged field of all the remaining electrons. The set of electron 
orbitals is obtained as the solution of this problem. Electrons are placed in this 
orbitals starting from the one with the lowest energy and filling the N/2 energetically 
lowest spatial orbitals where N denotes the number of electrons in the system. All 
the other orbitals remain unoccupied. The solution is obtained iteratively with the 
desired precision for the chosen basis sets. We do not want to discuss here the quality 
of results obtained with such a wavefunction as it is well known from both solid state 
physics and quantum chemistry that the Hartree-Fock approximation provides a 
reasonable but not accurate wavefunction and that for some properties it only gives 
qualitative results. However, it is also known that this wavefunction may serve as 
a good starting point for further improvement towards the exact wavefunction as 
long as weakly correlated systems are considered. So, our aim is start from the HF 
approximation and to go beyond it in order to obtain some good approximation to 
the correlated wavefunctions of a system and to calculate the ground state energy 
and energies of the excited states. 

Moreover, the system which we refer to is periodic and infinite. Our approach 
will allow to transfer the data obtained in clusters with standard quantum-chemical 
program packedges to the infinite system exploiting the local character of electron 
correlations. For this purpose we have to define local quantities (electron orbitals 
and matrix elements) which have identical meaning in the clusters and in the infinite 
systems and to esteblish a scheme to transfer them. 

In the first Section we consider the HF solution for infinite systems, namely, 
Bloch waves, the ground state energy and the band structure. In Section 2.2 we 
redefine these quantities in terms of Wannier orbitals. In Section 2.3 we show how 
those local quantities may be obtained for clusters and how they are transferred to 
infinite systems to give the band structure. In Section 2.4 we discuss why electron 
correlation should be taken into account, the physical meaning of this effect and how 
it is incorporated in our method. There we also define the correlated wavefunction. 
In Section 2.5 the quantum-chemical correlation method is described: why we have 
chosen the particular method, its advantages and disadvantages as compeared to 
other methods, the results which we get by using it and their quality. 
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2.1 The Hartree— Fock approximation 

The electronic ab initio Hamiltonian of a solid or molecule in the Born-Oppenheimer 
approximation is written in the form 

*..=-e^v?-j:-^+^^+^-^ (2.1) 

where r, and are the coordinates of electrons, and Za are the fixed coordi- 
nates and the atomic number of nucleus A. Atomic (Hartree) units are used here 
where electron charge e is set to 1. The last term in Eq. ()2.H) is the interaction 
energy between the fixed nuclei and does not affect the electronic eigenfunctions. 
We omit this term in further consideration bearing in mind that this constant for 
given geometry must be added to the system total energy at the end. 

The many-body wave function for independent electrons may be written as 
an antisymmetrized product (or the Slater determinant) of the spin orbitals ^/'^(ro"): 



$(ricri;...;r7vaAr) 



1 



V^i(rio-i) 
V^i(r20-2) 



^2(riO-i) 
V^i(r20-2) 



^l(r7va7v) ^i(rAraAr) ■■■ 

det[^/;i(riai) . . . 



V'Ar(riCri) 
V^Ar(r2Cr2) 

?/;Ar(riv(T7v) 



(2.2) 



Let us denote $(riO'i; . . . ; rTvcrjv) as $ for the convenience. The expectation value of 
the ground-state energy of the system is defined as Eq = ($|ifei|^)- By optimizing 
the spin orbitals ipi one obtains the lowest value for Eq under the constraint that spin 
orbitals ipi remain orthogonal to each other and normalized: (ipilipj) = 6ij. As the 
solution of this variational problem the HF single-determinant wavefunction of the 
many-electron system is obtained. It consists of spin orbitals being eigenfunctions 
of the self-consistent single-particle Fock operator F: 

Fm = e,\iJi) (2.3) 

where F takes a diagonal form in the representation of the canonical orbitals tpi and 
has the following matrix elements: 

N 

fij = + ^{{iJiiJilvuli^ji^i) - ii^ii^ilvuli^ii^j)) ■ (2.4) 

1=1 

and electron orbitals tpi are obtained self-consistently. In Eq. (j2.4|) index I counts 
the occupied orbitals and operators hi and f 12 are defined as 



hi 



Za 
r^i - ril 



(2.5) 
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and 

Vl2 = 7- (2.6) 

|ri - r2| 

The derivation of the equation ()2.3|) can be found in many textbooks (e.g., in |S()96j 1. 

Getting an analytic solution of the Fock equation ()2.3|) is not possible in general 
case and numerical methods have to be used for this purpose. In quantum-chemical 
program packages the Fock eigenfunctions are obtained iteratively in the form of 
linear combinations ipi = Xj'^ji within a chosen atomic basis sets {xj} by varying 
the coefficients aji and minimize the ground-state energy of the system. In our 
study solving the Fock equation is a routine procedure and we do not want to pay 
too much attention to the details of Hartree-Fock theory. We only want to point 
out the physical meaning of the eigenvalues in Eq. (j2.3j) . The quantity Si represents 
an energy of an electron occupying the spin orbital tpi. It consists of the electron 
kinetic energy, its potential energy in the field of nuclei and in the mean field of the 
other electrons. The latter in its turn consists of the Coulomb repulsion energy and 
the exchange energy which are determined correspondingly by the two last terms 
in the right-hand side of Eq. (j2.4j) . The spin orbitals with the lowest orbitals 
energies are occupied in the ground state. Further we use indices a, b, c, . . . for 
these occupied orbitals. The remaining eigenfunctions from ()2.3p stay unoccupied 
in the ground state and are labeled by r, s, t, . . .. If a statement or a formula is 
valid for both, occupied and virtual orbitals, we use the indices i, j, k, . . .. Also for 
the sake of convenience we use \i) instead of {ipi). Multiplying Eq. (j2.3|) by {i\ from 
the left and taking into account Eqs. (j2.4p - (j2.6p we get an expression for orbitals 
energies 



N N 

Si = {i\hi\i) + J2{{'b\vl2\^b) - {ib\vu\bi)) = + - {ib\M)) (2.7) 

6=1 6=1 

where the potential V12 is omitted in the matrix element for convenience. In partic- 
ular we can write the energies for the occupied orbitals as 

Ea = {a\h\a) + ^{{ab\ab) - {ab\ba)) (2.8) 
and that for the unoccupied orbitals as 

N 

= {r\h\r) + ^((r6|r6) - {rb\br)). (2.9) 
6=1 

These two quantities are of special interest in our study. Let us look more 
closely at expressions ()2.8p and ()2.9p . The energy of an electron in an occupied 
orbital ipa ()2.8|) consists of its kinetic energy, the energy of attraction to the nuclei 
plus the Coulomb {{rb\rb)) and exchange {—{rb\br)) interaction with each of the 
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N — 1 remaining electrons of the system in the ground state. Therefore we can refer 
to the value — as to an energy which is necessary to pay in order to remove an 
electron from the occupied orbital ipa without changing the rest of the system, that 
is the ionization potential of the system in the frozen orbital approximation. An 
unoccupied orbital energy is the sum of kinetic energy of an electron on this orbital, 
the energy of its attraction to the nuclei and energies of Coulomb and exchange 
interaction with all electrons of the system in the ground state as if an additional 
[N + l)-th electron would have been adiabatically added to the system (i.e. in the 
frozen orbital approximation) on the r-th orbital. Thus, we can regard the value 
—Er as the energy one gains when adding one electron to the neutral system and 
placing it into the orbital ipr that is by definition the electron affinity. These two 
statements are known as Koopmans' theorem |Ko34j . 

To demostrate how these simple relations ()2.8|) and ()2.9p arise from the frozen 
orbital approximation we provide the proof for this theorem. The ionization poten- 
tial is the difference between the energy of the system with — 1 electron (excited 
state) and the energy of A^-electron system (its ground state) 

IP= ^-^E-Eo (2.10) 

and the electron affinity is the difference between the ground state energy and the 
energy of (A^ + l)-electron system 

EA = Eo- (2.11) 

It is convenient to define state eigenf unctions in the second quantization formalism 
(especially in the case of infinite systems). Then the ground state wavefunction is 
written as 

\^) = ll^i\0) (2.12) 

a 

where c\ is the creation operator for an electron in an occupied spin orbital tpa and |0) 
is the vacuum state. The (A^ — l)-electron state in the frozen orbital approximation 
is given by 

I ^-1$,) = c,|$) (2.13) 

and the (A^ + l)-electron state by 

I ^+^,) = 41$) (2.14) 

where the operator Cc destroys an electron in occupied spin orbital ipc and the 
operator cj creates an electron in virtual spin orbital ipr- Let us note that ^~^c 
and ^~^^r are singe determinants constructed from A^ — 1 and A^ + 1 spin orbitals 
with 1 / y^(iV^^T)y and 1 / y^(iV-i-T)I as normalization factors, respectively. The 
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energies Eq, ^ ^Ec and '^'^^E^ are obtained as the expectation values of these single 
determinants: 

N N N 

Eo={<^\ H\ $) = ^{a\h\a) + -^^{{ah\ah) - {ah\ba)), (2.15) 

a=l a=l fe=l 

^-^E, = (^-^, I H I ^-^,) = Y.{a\h\a) + 1 ^^((a6|a6) - (a6|6a)) (2.16) 
and 



N 

N+^Er = \H\''+%.) = Y,{a\h\a) + {r\h\r) + 

a=l 

^ N N N 

2EE(<«^I«^) - {ab\ba)) +Y,{{rb\rb) - {rb\br)). (2.17) 

a=l b=l b=l 

The energy for taking an electron out from the spin orbital ipc of the neutral 
system is given by difference between (|2.16|) and (|2.15|) 

IP, = ^{a\h\a) + ^^^{{ab\ab) - {ab\ba)) 

N N N 

a=l a=l 6=1 

= -{c\h\c) - ^((c6|c6) - {cb\bc)) = -e, (2.18) 

and the energy obtained from adding an electron to the neutral system in the spin 
orbital ipr is given by difference between ()2.15|) and ()2.17p 

N 

EAr = -{r\h\r) - ^((r6|r6) - {rb\br)) = -e^. (2.19) 

b=i 

We have described the (A^ + 1)- and (A^ — l)-electron states by single-determinant 
wavefunctions constructed with frozen (non-relaxed) spin orbitals which are ob- 
tained as a solutions of the equation (j2.3j) for the ground state of electron sys- 
tem. For the case of finite molecules one should not expect spin orbitals of a neutral 
system to be the same as in cases of the charged one. For the more apropriate 
description of the (A^ + 1)- and (A^ — l)-electron states in molecules, the so-called 
ASCF approach, one has to solve the corresponding HF equations for the charged 
open-shell systems and get energies ^"^E, and ^^'^Ej. from these solutions. Thus, 
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the expressions derived above for IP,, and EAj. in the frozen orbital approximation 
only give approximate results for finite molecules. However, formulas ()2.18|) and 
(|2.19p are correct for extended systems where the electrons are delocalized over an 
infinite volume and hence no orbital relaxation takes place. Here we mean " correct" 
in the sense of the ASCF approach which is purely based on the HF energies. 

In a periodic infinite system the single-electron wave functions must obey the 
translation symmetry of the crystal. This leads to Bloch's theorem which states 
that each electron eigenfunction translated by a distance R = nai along a lattice 
vector should be equal to the original one within an exponential prefactor 

M^ + RW = e"'''Mr) (2.20) 

where, for simplicity, we only consider a one-dimensional crystal with lattice con- 
stant ai (aligned along the x-axis), n is an integer and the electron states are labeled 
by a crystal-momentum quantum numbers k G [—ir/afj+n/ai]. The fc-dependent 
eigenvalues of the states form energy bands and the quantity Sku refers to the en- 
ergy of an electron in a Bloch state with crystal-momentum k and from the band 
u. In the ground state of an insulator or a semiconductor the energetically lowest 
bands are completely filled (valence bands) and are separated by a gap from the 
others which are empty (conduction bands). For metals partially filled bands exist 
in the ground state. We explicitly exclude metals for our approach to correlated 
wavefunctions in solids and consider systems with band gaps only. 

As Koopmans' theorem holds for extended systems the electron energy Sku of a 
Bloch wave ipku is (up to a sign) the ionization potential for taking out an electron 
from the system (IPfci^) in case of valence electrons and the energy gain from putting 
an extra electron into a conduction band which is the electron affinity (EA^j,). As 
defined in (I2.1(J|) and (j2.1H) they are the differences of the energies of [N — 1)- 
or (A^ -|- l)-particle system respectively and the neutral A^-particle system. These 
energies are defined as the expectaction values of Hamiltonians for the corresponding 
single-determinant wavefunctions. Since A^, in this case, approaches infinity such a 
representation of the wavefunction becomes cumbersome and the use of the second 
quantization formalism seems to be more appropriate here. Then we define the 
ground-state wavefunction of the A^-particle system as 

occ 

i'^)=n^u^uio) (2-21) 

ku 

where operator c^^^ (c^j.) creates a spin-up (spin-down) electron in a Bloch wave, k 
runs over the first Brillouin zone and z/ counts only the valence bands. The energy 
of this state is Eq and is proportional to the volume of the extended system. The 
(A^ — l)-particle-state wavefunction is obtained by applying an annihilation operator 
for an electron with crystal-momentum k in the valence band i' and with spin a to 
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the wavef unction ()2.21|) 

\kiya) = Cku. I $) (2.22) 
and the (A^ + l)-particle-state wavefunction by applying a creation operator c^^^ to 

\kl^^)=cl,J^) (2.23) 

where the index is used to label conduction bands. In this new notation we can 
define bands as energy differences which depend on k. The valence bands are given 
by 

-IP;,,^ = -{kua I H I kua) + Eq (2.24) 
and the conduction bands by 

-EAfc^ = {k^ia I H I k^ia) - Eq. (2.25) 

At first glance, such a definition of the energy bands looks more complicated 
as compared to the energies of single-particle states in a crystal e^u and Sk^. How- 
ever, it will give us the possibility to express IP^u and EA^^ as a kind of Fourier 
transformation of some local matrix elements and to include electron correlation by 
affecting these matrix elements which would be not possible if electrons are treated 
as Bloch waves. 

Some alternative approaches to correlated band structure (e.g. based on DFT) 
are focused on finding an equation analogous to Eq. ()2.3j) but already containing 
ground-state correlations. In such schemes the band energies Skv are usually simply 
defined as the eigenvalues of a single-particle equation similar to the HF equa- 
tions. However, these values are not strictly energies of correlated holes or attached 
electrons and in the DFT framework this approximation usually leads to an under- 
estimation of the band gap of insulators and semiconductors. In contrast, in our 
method, we calculate separately and explicitly correlated (A^ — 1)-, (A^ + 1)- and 
A^-particle states and this allows us to see how the system reacts in detail on the 
presence of an extra charge. 



2.2 Wannier orbit als 

In the HF approximation an electron is considered as moving in the averaged self- 
consistent field of all the other electrons. But in reality each particular electron tries 
to avoid coming too close to its nearest neighbors due to Coulomb repulsion. This 
effect can not be described on the HF level and is called electron correlation. 

In our study we describe the process of adding of an electron to an infinite 
system or removing it which is equivalent to producing a positive charge. These 
extra charges are totally delocolized in the extended system since extra electrons (or 
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holes) in crystals are Bloch waves and the response of a system with infinite number 
of electrons to these extra particles is a pure correlation effect. Alternatively, we can 
use some mathematically constructed localized orbitals for the electrons and holes 
and restrict the response of an area around the localized extra charge (localized 
electron or hole). As a consequence, the correlation problem is reduced to the 
consideration of a limited number of electrons in a finite part of the crystal. The 
obtained correlation effect can subsequently be transferred back to Bloch states. 
This way we are able to get correlated wavefunctions and correlated band structures. 
In the present Section we will introduce the localized crystalline orbitals which can 
be used for that purpose. 

The natural choice would be the atomic basis functions however they are not 
mutually orthogonal in the crystal and are not solutions of the HF problem which 
we want to start from. The HF eigenfunction of the system, obtained in the form of 
the single determinant, in turn, has the nice property to be invariant (up to a phase 
factor) under any unitary transformation of the set of single-particle orbitals which 
compose it. Hence, a set of spatially localized orthonormal orbitals can be obtained 
as linear combination of canonical (delocalized) ones 



and several localization schemes (finding the unitary matrix uji which provides the 
less extended orbitals) have been proposed in quantum chemistry for finite molecules, 
e.g. |FB60j . |ER65j . |PM89j . Then, the wavefunction of the whole system can be 
written as a single determinant over (fi instead of ipi. The total energy of the 
system obtained as the expectation value of the Hamiltonian with the transformed 
wavefunction remains unchanged. 

In periodic systems the Wannier functions |W34j are known to yield localized 
electron orbitals if a suitable unitary transformation of the Bloch waves is used: 



and details of this step can be found in |E88j . The m Wannier functions (/?om which 
represent m energy bands z/i, . . . z/^ are referred to some unit cell (for instance by 
their centers) and repeated over the whole crystal by primitive translation vectors 
R. Pairs of spin-up and spin-down electrons from the m bands are said to occupy 
these repeated Wannier orbitals (WOs). Now the (A^ — 1)- and (A^-l-l)-particle-state 
wavefunctions are obtained from the canonical ones {kua) introduced in Section 2.1 
by a transformation totally analogous to the Wannier transformation of the Bloch 
orbitals. The wavefunction of the state with a localized hole is then given by 



(2.26) 



j 




(2.27) 



\Rna) = CRnal^) 



(2.28) 
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where an electron with spin a is removed from the n-th Wannier orbital centered at 
the unit cell with lattice vector R (we refer to a one- dimensional crystal here) and 
the ground-state wavefunction $ is constructed from Wannier orbitals. Note that 
due to translational invariance the vector R is not uniquely defined. However by 
fixing once the reference (or zero) unit cell we eliminate this ambiguity. Analogously, 
the (A^ -|- l)-electron states are defined as 

\Rma) = c^j,^J^) (2.29) 

where m is used to indicate the use of WOs from the conduction bands. 

Now we show how energy bands can be obtained from quantities defined in terms 
of \Rna) and \Rma). Analogously to our discourse in Section 2.1 we define "local" 
matrices of ionization potentials 

lPR>_ji^nn' = {Rna\H\R'n'a) - 6RR'6nnEo (2.30) 
and electron affinities 

EAR'_jirnm' = — {Rma\H\R! m' a) + 5RRi5mm' Eq. (2.31) 

Following the results (j2.18j) and (j2.19p their diagonal elements can be regarded as the 
energy of an electron occupying some Wannier orbital n (or m) in any cell and the 
off-diagonal elements are considered as hopping matrix elements between orbitals n 
in the cell R and n' in the cell R' . However, one should always keep in mind, that 
in reality an electron can not be removed from (or added to) some Wannier orbital. 
Nevertheless, these matrix elements can be used to obtain the band structure. 

Since the quantities ()2.30|) and ()2.3ip depend only on the difference R' — R we 
fix R at the zero cell {R = 0). In Fig. 12. II two occupied Wannier orbitals in trans- 
polyacetylene are shown schematically. The WO (pi is centered at the zero unit cell 
(marked by blue bonds). The energy of an electron in this orbital is — IPo,m. The 
hopping matrix element between WO yjj and WO ipj from the cell at distance R 
from the reference cell is — IP/j^jj. From now on we will refer to IP_R,ij and EA/j^jj as 
the "local matrix elements" (LMEs). 
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The energy bands (one-particle energies as a function of the crystal-momentum 
k) are obtained as eigenvalues of the "back-transformed" matrices 

lFnn'{k) = J2^"'''^^R,nn' (2.32) 
R 

and 

EA„„, (k) = J2 e'^^'EAR^mm' . (2.33) 

R 

The matrix elements lP/j^„„/ and EA/j mm' are the key quantities in our method 
for correlated band structures. Being calculated on the HF level they reproduce the 
HF band structure. In fact, only a relatively small amount of matrix elements is 
needed for doing this as the hopping matrix elements decay rapidly with increasing 
distance R and therefore the infinite summation in Eqs. (|2.32j) and (|2.33p can be 
terminated at some i?cut such that the rest of the sum is negligible with respect to 
the terms summed up so far. 

The effect of electron correlation leads to a change of the matrix elements IP R^nn' 
and EAR^mm' as will be shown in the following chapters. This means that the 
energy of a correlated hole or attached electron in some localized orbital and also 
the hopping matrix elements between different sites are not the same as in the HF 
approximation but contain all the information about the response of the charged 
[N — 1)- or [N + l)-particle system and the correlation in the ground state. Thus, 
substituting "correlated" matrix elements IPr^^^, and EA^^^, in Eqs. fl2.32|) and 
()2.33|) gives the band structure with electron correlations included. 

The correlation effect is determined by a variational method in which all unoc- 
cupied orbitals are used to construct the variational space. In a crystal there is an 
infinite amount of virtual orbitals which would lead to an infinitely large variational 
space and make the calculation not feasible. To arrive at a finite variational space 
all the unoccupied Bloch waves can be transformed to Wannier orbitals and ordered 
by their distance from the localized charge and then only WOs up to a certain 
distance -Rcut are taken into account. However, firstly, the localization (Wannier 
transformation) of all virtual Bloch orbitals becomes very expensive for rich basis 
sets and, secondly, WOs corresponding to energetically high Bloch states are spa- 
tially rather extended and therefore the cut-off radius i?cut has to be chosen large 
enough that may again lead to infeasible correlation calculations. As an alternative 
to the set of unoccupied WOs, one can use atomic orbitals with one projected out of 
the occupied space. This idea is being efficiently implemented in quantum-chemical 
codes (see e.g. |SHW99j and |WMK03] ). These orbitals are perfectly localized and 
the cut-off radius is reduced to the minimum value but they are (i) not mutually 
orthogonal and (ii) can not be used to reproduce the HF bands by a suitable unitary 
transformation. 

In our approach we avoid the problem of an infinite variational space by per- 
forming the correlation calculations in finite molecules representing a small piece of 
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the periodic system. This imphes a finite number of virtual HF orbitals in these 
molecules. 

2.3 One-particle configurations 

Here we have come to the crucial point of our method. The matrix elements IPr^ij 
and EAij jj defined in fl2.3()j) and 1)2.311) can be obtained (approximately) from calcu- 
lations in finite molecules. Moreover, electron correlation can be easily incorporated 
by performing standard quantum-chemical finite-cluster correlation calculations. In 
the present Section we focus our attention on the first statement. 

For understanding this statement the molecular-orbital theory may provide a 
hint. A bond orbital (fa between two atoms A and B in a. molecule (a localized one- 
particle wavefunction) is obtained as a linear combination of atomic basis functions 
with the main contribution centered on these two atoms: 

Va = '^ aA,iXA,i + ^ 0!B,iXB,i + ^ ac,iXc,i (2.34) 

i i i,C^A,B 

with |ayi,i|, |a_B,.j| ^ c^a b l^c,j| where XA,i is the i-th atomic basis function of atom 
A and \ac,i\ decays rapidly with the distance between the C-th atom and the bond. 
The same bond one expects to find in a covalent extended system when the same two 
atoms at the same distance and with similar atomic structure in the neighborhood 
are considered. Therefore one expects to obtain occupied Wannier orbitals in a 
crystal which are very similar to the localized occupied spin orbitals in a molecule 
when this molecule represents some finite part of the crystal provided the same 
atomic basis sets and the same localization criteria are used for both the molecule 
and the crystal. The larger the molecule the closer the localized molecular orbital 
(LMO) in the central part of the molecule are to the relevant WO of the crystal. 
The LMOs at the edges of the molecule are less similar to the corresponding WOs as 
the surrounding of the former differs strongly from the surrounding in the crystal. 

The matrix elements IPr^uu' ()2.30p are defined in terms of localized crystalline 
orbitals which are derived from HF solutions for the closed-shell system in the ground 
state. The diagonal elements — IPo,nn are the energies eo.nn of an electron in the n-th 
orbital of the reference unit cell which constitutes a generalization of Koopmans' 
theorem to the case of localized orbitals. The quantity £o,nn depends on the shape of 
the orbital itself and on the atomic and electronic structure of the close surrounding. 
If they are similar in the molecule and in the crystal the corresponding orbital 
energies should coincide. The same statement is valid for the hopping matrix elements 
in ()2.3(Jj) when the two localized orbitals and the structure between them do not 
differ. Thus all the matrix elements which are needed in the restricted summation 
()2.32|) can be obtained from calculations in large enough but finite molecules. 

Similar argumentation holds for the unoccupied WOs corresponding to the en- 
ergetically lowest conduction bands with the distinction that localized unoccupied 
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orbitals in molecules can not be obtained by standard quantum-chemical localiza- 
tion methods, in general. However, in this case the WOs from the periodic system 
can directly be projected onto a molecule when the structure and the basis sets are 
identical and can be used as unoccupied LMOs. This technics is described in details 
in chapter 3.2.2. Here we only emphasize that local matrix elements EAn^rnm' ()2.3H) 
can be obtained from calculations in finite molecules as well. 

To get the matrix elements ()2.30p and ()2.31|) from molecular calculations an 
appropriate molecules have to be chosen. The structure of such molecules must 
contain the appropriate clusters of a crystal under consideration, e.g. those unit 
cells which are needed to represent the corresponding WOs without noticeable loss of 
accuracy. Such clusters cut from covalent crystals have dangling bonds which can be 
saturated by hydrogen atoms. To be more specific let us consider the occupied WOs 
schematically shown in the Fig. 12. II There, four unit cells of trans-polyacetylene are 
drawn. The two WOs define the three matrix elements IPo,m, IP_R,ij and IPnjj = 
IPojj. To obtain these matrix elements from a molecule calculation the shown 
cluster CgHg (with four unit cells) is appropriate and therefore for calculations one 
can use a CsHio molecule with the same geometry as the CsHg cluster but with 
the two cut C-C bonds being substituted by C-H bonds with the same orientation 
as the C-C bond and a typical C-H bond length. In our investigation of trans- 
polyacetylene we use CgHg, CgHio, C10H12 and C12H14 molecules which represent 
different-size clusters of the infinite chain all being constructed to contain full unit 
cells such that they can be terminated by single C-H bonds. The size of the chosen 
cluster is defined by the distance between the relevant WOs and should be as small 
as possible for cheaper calculations but still big enough to preserve the shape of the 
WOs and by this the accuracy of the results. These criteria hold for the selection of 
proper molecules for all kinds of covalent materials, e.g. diamond, silicon, gallium 
arsenide, polyethylene etc. 

Thus, having chosen appropriate molecules we can formally redefine the matrix 
elements (j2.30|) and (j2.31|) in terms of localized molecular orbitals 



IP 



mol 
IJ 



imol 




(2.35) 



with 



I J) = CjI^'""') 



(2.36) 



and 



EA 



mol 

■I*J* — 



{riH'^^y*) + Si,j,Ei 



rmol 



(2.37) 



with 



r) = 4.1$"^°^) 



(2.38) 
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where c/ destroys an electron in an occupied localized molecular orbital named I (a 
bonding orbital) and c|, creates an electron in an unoccupied localized molecular 
orbital /* (an anti-bonding orbital). 

Starting from this point we shall refer to local matrix elements as obtained from 
calculations in molecules. Also we adopt the notations for many-body wavefunctions 
used in quantum chemistry. In the HF approximation the electronic wavefunction 
of a molecule in the ground state is written as a Slater determinant ()2.2|) . Since the 
Slater determinant is essentially (up to a phase factor) invariant under any unitary 
transformation of the orbitals we rewrite it in terms of the localized spin orbitals 

<^ = -^det[ipu....^N] (2.39) 

where is the number of electrons in the molecule. In the HF ground state config- 
uration all occupied spin orbitals (pa contain an electron and all virtual spin orbitals 
iPr are empty. By taking out an electron from an orbital (fa we produce an (A^ — 1)- 
particle state with the wavefunction 

= , det[(y9i, . . . , Lpa-1, <fa+l, • • • , ^n]- (2.40) 

a/ (A — 1)! 

When an electron is added to the system and placed into the virtual spin orbital (pr 
an (A^ + l)-particle state is produced and its wavefunction is 

'^r = ^=^=det[<pi,...,^N,^r]- (2.41) 

^(iV + Ij! 

We shall call the wavefunctions ()2.40|) and ()2.41|) "one-particle configurations" as 
they correspond to states with one particle (a hole or an electron) being added to 
the system in the ground state. 

Now we would like to show explicitly how local matrix elements are obtained 
from molecular quantum-chemical calculations on the HF level. As follows from 
Eqs. ()2.35p - ()2.4H) the LMEs are given in terms of the one-particle configurations as 

IPab = ($.|i/°^°l|<l>6) - SabEr\ (2.42) 

EA,,, = -(^.lif'^^'l^,) + 6rsEr' (2.43) 

where 

= (<l>|i/"^°i|<l>). (2.44) 

In the following we omit the superscript "mol" since our calculations are always 
performed in molecules except of the localization of the unoccupied orbitals and 
there we explicitly distinguish infinite and finite systems. From the HF calculations 
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with the program package MOLPRO |M2nnnj the total ground-state energy Eq, 
the canonical HF spin orbitals ipi and the Fock matrix ()2.4j) are obtained. In the 
representation of the canonical orbitals the Fock matrix has diagonal form 

= {^,\F\^,) = 6,,e, (2.45) 

where are the canonical orbital energy (j2.7|) . By a localization procedure applied 
to the canonical orbitals one obtains the LMOs and the corresponding unitary 

matrix Uji. Then using Eqs. fj2.26|) . ()2.16|) . ()2.10|) and ()2.18p one can express matrix 
elements IPab in terms of calculated quantities: 



occ 



i j 
occ occ occ occ 

i 3 i j 

occ occ 

= y^^UigUibjEo - Ej) - 5abEo = UiaEjUib. (2.46) 

i i 

(Here we omitted complex conjugation as molecular HF orbitals are usually real.) 
Thus, the local matrix elements IPab obtained on the HF level are in fact the elements 
of the Fock matrix (with the opposite sign) which is defined in terms of the localized 
occupied orbitals: 

lPab = -fab = -{Va\F\Vb). (2.47) 

Similarly one gets 



vir 



EA,.^ = - ^ UirUis{Eo + Ei) + SrsEo = - ^ UirEiUis- (2.48) 
i i 

which simply is 

EA,, = -frs = -{Vr\F\^s). (2.49) 

In ()2.46|) index i runs over all occupied canonical orbitals and in ()2.48|) it counts 
only relevant virtual ones which means all virtual canonical orbitals which were used 
for the generating of the virtual LMOs ipr- 

The crystalline matrix elements (j2.3(Jj) and (j2.3ip which enter the restricted sum- 
mations (j2.32p and (j2.33|) can be replaced by the matrix elements (j2.46p and (j2.48j) 
obtained in molecules, and by this the HF band structure of the corresponding crys- 
tal can be calculated. Such a band structure, obtained by LMEs from molecules, 
can be compared to that obtained from canonical HF calculations for the infinite 
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system, e.g. done by the program package CRYSTAL |CRY98j using identical ba- 
sis sets and geometries of the crystal and the clusters. Also the LMEs themselves 
can be compared directly if elements of the Fock matrix in the representation of the 
WOs can be obtained from the calculations of the infinite system. Such comparisons 
show how accurate the cluster approach is. In fact, there should be small devia- 
tions between the band structure by LMEs on the HF level and that calculated by 
CRYSTAL because of the truncation in the summations ()2.32|) and ()2.33|) . When 
the correlated band structure is obtained by LMEs including correlation corrections 
one can account for these small deviations as is explained in chapter 4.2. 

2.4 Electron correlation effects 

Above we have showed how the many-body wavefunction in the HF approximation 
can be obtained. This approximation serves as a good starting point for the study 
of electronic properties both in finite molecules and in extended systems providing 
single-electron orbitals with corresponding orbital energies, the ground state energy 
of the system, the energies of the {N + 1)- and {N — l)-particle system in different 
excited states, ionization potentials and electron affinities. However, all these results 
can only be considered as approximate since the effect of electron correlation is 
missing on the HF level. 

Let us consider an extended periodic system (a crystal or a polymer chain) in 
the ground state for which the HF solution is known, in particular the ground-state 
energy Eq and the localized electron orbitals (Wannier orbitals). Eq accounts for 
the Coulomb repulsion between the electrons only in an averaged way (mean-field 
approximation). In reality electrons move in accordance to the movement of their 
neighbor electrons trying to avoid coming too close to each other and to reduce by 
this the Coulomb repulsion. Consequently, the total energy of the system, Eq, is 
reduced noticeably. Alternatively, we can describe this effect as charge fluctuations. 
In a neutral system a dipole may emerge due to local fluctuation of the charge. 
The system reacts to this by inducing another dipole and gains the energy of the 
dipole-dipole interaction. 

If an additional electron is added to a localized virtual orbital (or removed from 
some occupied WO) the system reacts on the presence of the extra charge by pro- 
ducing a polarization cloud around it. The main effect of electron rearrangement 
happens in the close surrounding and decays with increasing distance. Of course, 
the effect, described above for the neutral system, also takes place in the charged 
one. Both these effects lead to a substantial reduction of the energy ^~^^E of the 
(A^ -|- l)-particle system and the energy ^~^E of the (A^ — l)-particle system. 

As we have seen, electron correlations in the ground state of a system are of 
van der Waals type and the contribution to the system energy falls rapidly with the 
distance between dipoles. Therefore, only rather small clusters are needed to prop- 
erly estimate the correlation contribution to the ground-state energy. Signiflcantly 
more efforts requires the proper description of electron correlations in the charged 
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system. An accurate treatment of the polarization cloud around the localized extra 
charge is necessary as this is the main contribution to the correlation corrections to 
the LMEs. Since electron-dipole interaction decays more slowly than dipole-dipole 
interaction, larger clusters have to be chosen. Also, one can account separately for 
the (relatively small) contribution coming from the polarization of the farther sur- 
rounding which is not present in the cluster models. For that purpose the crystal 
is considered as a polarizable continuum in the presence of a point charge and the 
approximate formula for the corresponding contribution to the correlation energy 
is presented in chapter 4.3. Thus, electron correlations have mainly local character 
and the corresponding corrections to the ground-state energy and the LMEs can be 
obtained from calculations on finite clusters. 

Now we would like to show explicitly how electron correlation is included in 
quantum-chemical post-HF methods. Solving the Hartree-Fock equations ()2.3|) - 
fl2.(j|l not only provides an approximate ground-state wavefunction in the form of 
a single determinant ()2.2j) but also one obtains a set of 2 x A^bas spin orbitals (or 
A'bas spatial orbitals) where A'bas is the number of atomic basis functions {Xi}, 
i = 1, . . . , A^bas for a given molecule. Using these spin orbitals one can select a 
huge number of different orbital configurations which is of the order of (^''''°) where 
Q) is the binomial coefficient. Each of these orbital configurations represents a 
single determinant referred to simply as " configuration" . The HF wavefunction $ is 
usually defined as that configuration which gives the lowest (self-consistent) energy 
Eq for a given number of electrons in the molecule. 

To incorporate correlations in the wavefunction one takes advantage of the fact 
that the entire set of configurations $/ form a basis of the Hilbert space. Hence one 
can write the "correlated" wave function of the system as the linear combination of 
single-determinant configurations: 

$-- = $ + ^a^$^ (2.50) 

(here written in the intermediate normalization without a prefactor for $). By vary- 
ing the coefficients «/ one can minimize the energy corresponding to the correlated 
wavefunction (|2.5Up 



^^corr 




|/,|,corr , 





(2.51) 



Of course, in practice, one has to truncate the set of configurations somehow. This is 
the essence of the so-called "configuration interaction" (CI) methods. The obtained 
energy Eq°^^ is lower than Eq from the HF calculations. The same argumentation 
holds for the wavefunction and the energy of (A^ — 1)- and (A^ + l)-particle systems 
when the (A^ ^1)- and (A^ + l)-electron configurations are used, respectively, to 
construct the multi-determinant wavef unctions. 

In fact only for small molecules with a few atoms and small basis sets a variation 
of all configurations is possible (full CI). For larger systems a drastic reduction of the 
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variational space (the number of involved configurations) must be done. The excited 
configurations denoted as $/ in ()2.50p can be sorted by the number of substitutions 
of occupied spin orbitals in the ground-state configuration by unoccupied ones and 
the correlated wavefunction ()2.5()|1 can be rewritten as 

*^°" = * + E«+ E E <tcKfc + --- (2.52) 

ar a<b, r<s a<b<c, r<s<t 

where and are the excited configurations, subscripts denote the spin 

orbitals excluded from $ and superscripts denote the substituting virtual spin or- 
bitals. 

The contributions to the correlated wavefunction, associated with configurations 
with three and more substitutions (so-called excitations), are negligibly small all 
together compared to those of first three terms on the right-hand side of Eq. ()2.52|) . 
Therefore, we take into account only the ground-state configuration (the dominant 
contribution to the correlated wavefunction) and the singly and doubly excited con- 
figurations (single- double CI). Note that the singly excited configurations do not 
contribute directly to the correlation energy of the ground state as they do not cou- 
ple to the HF ground-state configuration: = 0. This is the statement of 
Brillouin's theorem (see, e.g., chapter 3.3.2 in |S096j ) which holds as long as the HF 
spin orbitals are determined self-consistently. Thus, the second sum in Eq. ()2.52p 
including the double excitations gives the main contribution to correlation energy 
and includes dipole-dipole interactions. 

The correlated wavefunctions for the charged systems can be constructed simi- 
larly. For instance, the wavefunction of a system with a static hole in a frozen orbital 
ipa would be written as 

*r = <f a + 5^ + ^c^a^'c (2.53) 

b^a,r b<c^a, r<s 

Here both sums give contribution to the energy of the (A^ — l)-electron system. Ma- 
trix elements ($^^|if|$a) do not vanish and describe the polarization cloud around 
an extra charge. 

Of course, in reality a hole is delocalized over the entire system and a more 
general ansatz for the correlated wavefunction of the (A^ — l)-particle system has to 
be made: 

N-l^corr ^ ^ ^^^^ ^ ^ ^j^^j^ ^ ^ ^j^s^^rs (2.54) 

a a,b^a,r a,b<cj^a, r<s 

The number of configurations (the number of single determinants) required for 
this is huge even though we restricted ourselves to singly and doubly excited configu- 
rations only. However, we do not need to store the explicitly correlated wavefunction 
consisting of billions of single determinants since we are rather interested in physical 



26 



CHAPTER 2. LOCAL APPROACH TO WAVEFUNCTIONS 



quantities which are obtained with this wavefunction, namely, LMEs hke the ones 
introduced on the HF level ()2.42j) and ()2.43p . Below we show how this is done. 

First, we turn back to the representation of the HF one-particle configurations 
in terms of canonical orbitals in which the Fock matrix is diagonal. We denote 
these configurations as \E'j here with an electron being removed from (or added to) 
the canonical spin orbital ipi. The configurations \E'j can also be obtained as the 
linear combinations of the "local" configurations $j: \E'j = Ylj^j because 
the canonical spin orbitals can be obtained from the localized ones with the help of 
the unitary matrix Uji from the localization procedure (j2.26|) : ipi = Ylj'^j 
The correlated counterpart of the HF wavefunction \E^j is that eigenstate of 
the charged system which resembles the HF configuration \Efj the most. It will 
essentially be composed by local one-particle configurations like the ones in the first 
sum on the right-hand side of Eq. ()2.54|1 for ^-i<|)Corr rjj^g wavefunctions \E'™" of 
the {N — l)-electron or (A^ -|- l)-electron system are associated with total energies 

J^COTT . 

^^corr _ ^corr^corr_ ^2.55) 

Next, we formally introduce a one-to-one correspondence between the HF many- 
particle wavefunctions and their correlated counterparts by the use of a wave oper- 
ator Q which produces the correlated wavefunctions vl/™'"'' from the HF ones 

= Q^i, i = l,...,d. (2.56) 

Then, we divide the full space of local configurations ^j, and ^j^l used to 
construct the correlated wavefunctions vE'™'''' into a small subspace 7 of dimension 
d consisting of all SCF one-particle configurations (or equivalently of all local 
configurations $j) with P being the projector onto this space and its orthogonal 
complement Q (with the projector Q), so P + Q = 1. The J-space is often called 
the model space (or the space of model configurations) and the Q-space is called the 
external space, and accordingly the configurations which belong to the Q-space are 
called external configurations. Thus, the wave operator Q acts on configurations 
\E'j from the CP-space and yields the exact wavefunctions from the full Hilbert 
space. 

The HF wavefunctions \E'j can be obtained from the correlated counterparts by 
the back transformation 

^i = U (2.57) 

H is an operator acting from the dimensional space spanned by the correlated 
counterparts {'^i°"}i=i^...,d to the model space CP. It is similar but not identical to 
the projector P which can be written as 



d 

P = 5^ I VI/,) (VI/, 

1=1 



(2.58) 
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Figure 2.2: Schematic representation of the correlated wavefunction. 



because the projections P'^^°^^ are neither perfectly normalized nor mutually or- 
thogonal. In Fig. 12.21 we schematically show a correlated wavefunction Both, 
its projection PvE'™'''' and the HF counterpart lie in the ^-space but differ slightly. 

Let us substitute (j2.56j) into the left-hand side of Eq. (j2.55p and act on the 
obtained equation by the operator 11 from the left 



UHQ'^i = ni/^J°" = E^°''U'$'r^'' = (2.59) 

This way we get an eigenvalue equation for and E^°". Inspired by this relation 
we define an effective Hamiltonian as 



H"^ = UHn. (2.60) 
It operates fully in the IP-space and its eigenvalues are those of the exact Hamiltonian 



if^ff^. = E^°''-(I/i. (2.61) 

As the unitary transformation ()2.26p . which gives the localized HF spin orbitals 
of a molecule in terms of the canonical ones, is known from the localization procedure 
we can readily get the matrix elements of if'^^ in terms of the localized one-particle 
configurations 



i' j' 
= ^i'J- ^^'i' = Yl ^^''^ (2.62) 



where the energies E^,°" are those obtained in our correlation calculations. To be 
more precise, these energies as well as the elements of the matrix Uij were calculated 
by the MOLPRO program package |M2nnnj . |WK88j . |KW88j . |KW92j . The matrix 



elements Hfj^ defined in Eq. fj2.62j) differ from the corresponding ones on the HF 
level Hij and contain electron correlation as far as it is included in the energies 
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^corr_ Together with the correlated energy of the ground state Eq°" ()2.51|) the 
matrix elements of the effective Hamiltonian define the correlated LMEs: 

IPT = H^f - SabEr (2.63) 

and 

EA^r = -H^f + SrsE'o°'' (2.64) 

which are to be used in ()2.32j) and ()2.33|1 instead of the SCF ones to give the 
correlated band structure. Note the similarity of the equations ()2.fi2|l - ()2.64|) with 
the corresponding ones for the matrix elements IPa^ and EA^^ on the HF level ()2.46p 
and (E3H1)- 

It should be mentioned here that one can also use P instead of 11 and a corre- 
sponding wave operator W (which is linked to Q hj W = Q{PQ)~^) to define an 
effective Hamiltonian 

H"^ = PHW. (2.65) 

Despite the fact that H^^ is non-Hermitian it has the same properties as H^^ . 

At the end of this section we would like to say a few more words about the wave 
operator Q and its inverse 11. They enter the definition of the effective Hamiltonian 
()2.60|) however they are not used for obtaining its matrix elements H°^. The latter 
ones are calculated using equation ()2.62j) . Hence, Q and H (or P and W) are only 
introduced formally in our method. Nevertheless, P and W can be constructed 
by the use of perturbation theory up to any desired accuracy as is comprehensively 
explained in chapter 9 of Ref . |LM86j . Then the correlated wavefunctions ()2.54|) can 
be obtained by applying W to the projection Pvl/™'''' which is obtained perturbatively 
as well. 

2.5 Multireference configuration interaction method 

Above we have used the CI ansatz to write the correlated wavefunction of a many- 
body system ()2.52|) and ()2.54|) . However, there are also other correlation methods 
which use the HF wavefunction as a starting point and excited configurations to 
construct its correlated counterpart. Here we would like to mention two of them, 
namely the single and double excitation coupled-cluster (CCSD) method and the 
second-order M0ller-Plesset perturbation (MP2) method. The CCSD method is 
more successful in ground-state correlation calculations than the single and double 
configuration interaction (CI(SD)) and both have the advantage that they are size- 
consistent in contrast to CI(SD). Unfortunately, CCSD and MP2 are not suitable 
for excited-state calculations. 
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The CCSD ansatz for the correlated ground-state wavefunction is 

<^'Sh' = exp(j2t:Cacl+ t:iCaC,clA<^ (2.66) 

^ ar a<b, r<s ^ 

where the operators Ca and Cb (cj and cj) destroy (create) an electron in an occupied 
(unoccupied) spin orbitals v?a and (pb {^Pr and ips) of the ground-state configuration, 
respectively. The coefficients and t^l in Eq. ()2.66p are obtained from a set of 
algebraic equations. The exponential form of ^c!C iniplies that the correlated wave- 
function of a system consisting of some well separated subsystems is the product 
of the wavefunctions of these subsystems which is the crucial prerequisite for size- 
consistency. There is no simple extension of the CC method to multireference (MR) 
calculations which are needed to describe excited states. Though such an extension 
of CC became a hot topic in the literature in the last decade, e.g. |PB99j . |AMI00j . 
P03j . the MRCC method seems still to be far from the stage of implementation in 
standard quantum-chemical program packages at the present moment. 

The second-order M0ller-Plesset perturbation method can be generalized to the 
multireference case. However, in the perturbation expansion for the correlated wave- 
functions and energies of the {N — 1)- and (A^ + l)-electron systems there are terms 
which are proportional to 

(2.67) 

TP TP r ^ ' 

Here Ei is the SCF energy of some hole (or attached electron) state and E-^ is the 
energy of an SCF configuration with one hole (or attached electron) state plus one 
electron-hole excitation. It may happen, now, that E^ gets very close to E^^ in the 
case where E^ corresponds to an energetically deep lying hole and E-^ corresponds to 
an energetically high lying hole and a close electron-hole excitation. This situation 
takes place in systems with relatively small band gap and a strong dispersion of the 
bands. Such small denominators blow up the correlated energies and lead to wrong 
results. Therefore, the MP2 method can only be applied to systems with big band 
gaps and rather fiat bands |ARM98j . Also deep valence bands and high conduction 
bands can not be treated properly. This problem was discussed in, e.g. |SB96j and 
lAKSOlj where the MP2 method was generalized for crystal-momentum space and 
could be applied only to the vr bands of irans-polyacetylene single chains. 

As a variant, a hybrid method was used to determine the correlated valence 
bands of diamond, silicon and germanium in Ref. inEEna, f^m7\ and jXFMI] 
where single-excited configurations were incorporated by MRCI(S) and for the dou- 
ble excitations a separate calculation was performed employing quasi-degenerate 
variational perturbation theory |CD88j . In that case small denominators do not 
show up. However, the size-consistency problem of MRCI(S) is still present (though 
reduced as compared to MRCI(SD)) and the contribution of coupling terms between 
singly- and doubly-excited configurations to the correlation energy is missing totally. 
Also such a hybrid method is technically more complicated than pure MRCI(SD). 
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Hence, for our correlation calculations we have chosen the multireference version 
of the single and double configurations interaction method (MRCI(SD)) to calculate 
the correlated ground-state energy (|2.51|) and the matrix elements of the effective 
Hamiltonian ()2.(i2j) on both HF and correlation level. The size-consistency problem 
has been accounted for by the development of a special size-consistency correction 
(for details see Section 3.5). 

Let us first consider the CI(SD) ansatz for correlated ground-state wavefunction 
^corr_ r^YiQ HF ground-state wavefunction $ is written as the Slater determinant 
composed from localized occupied spin orbitals. In general all singly- and doubly- 
excited configurations are provided to the CI(SD) expansion. However, one can 
further restrict the number of excited configurations by forbidding to take out an 
electron from some particular spin orbitals. For instance, we can "freeze" the arti- 
ficial C-H bonds which are used to saturate the dangling bonds of the clusters. By 
this technique we exclude artificial contributions to the correlation energy associ- 
ated with these C-H bonds. Also, we can freeze the majority of the bond orbitals 
of a molecule and allow excitations from a few particular bonds only. By this, we 
can replace an expensive correlation calculation for a big molecule by a sequence 
of cheaper ones and sum up, in a proper way, the correlation contribution to the 
ground-state energy coming from the several parts of the molecule. This is the idea 
of the method of local increments introduced in |St92aj and |St92bj which we have 



adopted and which will be discussed in more details in Section 3.4. 

Thus, the correlated ground-state wavefunction of a molecule with electrons 
and n bonds (or 2n spin orbitals) open for constructing excited configurations is 
written as 



2n 2n 



= $ + ^ ^ a:<i>: ^I'^ai- (2.68) 



a=l r a, 6=1 r,: 



The coefficients and are determined by minimization of the ground-state 
energy of the system given by Eq. ()2.5H) where they play a role of variational pa- 
rameters. The resulting eigenvalue problem can be obtained by multiplying the 
Schrodinger equation 

^|$corr^ = ^corrj^corr^ ^2.69) 

successively by ($|, ($*| and ($*^| from the left: 



a,b,r,s 

i \ -r -r I T T^ it \ r 



a,r a,b,r,s 



, rs 



7-icorr 
= ^0 


(2 


70) 


T-icorr^t 
- ^0 «c 


(2 


71) 


TTicorr „ tu 
= ^0 ^cd- 


(2 


72) 
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This eigenvalue problem is too large to be solved directly and an iterative procedure 
(according to Davidson |D75p is employed. The correlated ground-state energy 
can be obtained with any desired precision in this procedure. In our calculations 
the threshold for energy change form one iteration to the next is set to 10~^ eV. 
However, this does not mean that we get the full correlation energy with such an 
accuracy. That would only be the case for a full CI calculation (with all possible 
excited configurations included). Rather it defines an upper bound for the exact 
correlation energy since truncated CI is a variational method. 

Now we would like to show how the CI(SD) ansatz can be extended to get the 
matrix elements of the effective Hamiltonian ()2.55|) . Suppose we want to get energies 
of n states with an electron removed from one of the n bonds of a molecule and also 
the hopping matrix elements between these bonds. As bonds are spatial orbitals 
containing two electrons with opposite spins per bond in the ground state, two 
spin-degenerate states are possible for a hole in some particular bond. To resolve 
the problem we consider only states of a particular spin symmetry, e.g. only one 
particle configurations obtained by annihilation of a spin-up electron. Then we have 
to provide n corresponding spin-adapted configurations $i {i =1, 2, . . . , n) into the 
model space. They define an n x n matrix Hij = on the HF level. The 

model wavefunctions $j are first transformed to diagonalize this matrix: 

n n 

The configurations '^i are called reference configurations or reference states and Ei 
are the energies of these states. In our calculations we get the SCF energies Ei and 
coefficients Uij from a preliminary step of the MRCI calculations. Then the SCF 
matrix elements are obtained as 



Hij = ^UkiUkjEk (2.74) 

k=l 

in close analogy to Eq. ()2.62p for -fff/ • 

To produce excited configurations for the CI expansion of the correlated wave- 
functions ^E'i"", pair excitation operators CacJ and CaCbclcl are applied to each of the 
reference states \E'j, yielding the following ansatz for the n correlated wavefunctions: 

*r = E f^^^^) + + T!^/c.mcaCbcic^)j (2.75) 

j=l ^ a,r a,b,r,s ' 

which (in case of an (A^ — l)-particle system) is equivalent to 

= E {^^^)^^ + Y^<^)^^a + Y^o^-im-ai) (2.76) 

j=\ ^ a,r a,b,r,s ' 
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since the n reference states and the n spin-adapted configurations $j span the 
same (reference) space. In Eq. fl2.76p prime over sums means that electrons are 
annihilated from the n active bonds and only the remaining spin-down electron can 
be removed from the localized orbital j. These n eigenfunctions ()2.7(j|l correspond 
to n eigenvalues of the full Hamiltonian: 



and djabi''') ^i-re determined iteratively from a set of equations obtained by multiply- 
ing both sides of the n equations ()2.77|) successively by and from 
the left. For details we refer to |iWK88) . 

Once the iterations are converged the n correlated energies of (A^ — l)-electron 
system Ei°" can be evaluated. Being used in ()2.74|1 instead of the HF energies Ek 
they determine the matrix elements of the effective Hamiltonian H^^ as have been 
discussed in the previous section (see Eq. ()2.62p ). 

In the case of attached-electron states the wavefunctions \E'j in Eqs. (j2.75p corre- 
spond to the m (A^+l)-particle states given by linear combinations ()2.73|) of the local 
(A^ + l)-particle configurations with an extra (spin-down) electron being added 
to one of the m localized virtual orbitals. The MRCI-ansatz for the m correlated 
wavefunctions is written in a form similar to Eq. (|2.76|) : 



and here primes mean that only spin-up electron can be put into a localized virtual 
orbital (fi if it is already occupied by the attached spin-down electron. The obtained 
m correlated energies of the (A^ + l)-electron system El°" can be used in Eq. ()2.74j) 
instead of the corresponding m HF energies to calculate the matrix elements of the 
effective Hamiltonian Hf^ for the case of attached-electron states. 

Thus, equations (j2.73|) - (j2.78p define how the MRCI(SD) method is used to eval- 
uate the matrix elements of the effective Hamiltonian in the case of both (A^ ~ 1)- 
and (A^ + l)-particle states. 





m. 



(2.78) 



Chapter 3 

Realization of the method 



While we have outhned the general idea of our approach to the correlated band 
structures in the previous chapter we describe in detail the actual realization of this 
method here. The concepts we have developed to make such calculations possible 
are presented in the given chapter. We will refer to trans-polyacetylene (tPA) single 
chains as a test system. In section 3.1 the geometry and electronic properties of 
tPA on the HF level are described. In section 3.2 we show how to generate suitable 
localized occupied and virtual HF orbitals in the clusters. The corresponding HF 
local matrix elements are used to calculate the SCF band structure of tPA. In 
section 3.3 we discuss how to evaluate the correlation corrections to the LMEs by 
the MRCI(SD) method and the use of an approximate but well controlled scheme 
for this is presented in section 3.4. The possibility to reduce the size-extensivity 
error of the MRCI(SD) method is explored in section 3.5 and the analytic formula 
for the corresponding correction to the MRCI(SD) correlation energy is derived for 
open-shell systems. The occurrence of satellite states in cluster models and the 
handling of this problem by an appropriate construction of the model space in the 
MRCI calculation are discussed in section 3.6. 



3.1 The investigated system 

In Chapter 2 we have described the general idea of our wavefunction-based method 
for obtaining band structures of insulators and semiconductors with electron corre- 
lation effects being properly included. This approach was used in Ref. |GSF97j and 
|AFS00j to get the correlated valence bands of three structurally similar semicon- 
ductors (diamond, silicon and germanium). In the present work we make the first 
attempt to treat systematically not only the valence but also the conduction bands 
of a periodic system. For this purpose we introduce localized virtual HF orbitals 
which are, in contrast to localized occupied orbitals, usually extended over several 
unit cells of the crystal. Therefore, rather big clusters have to be used for a proper 
treatment of the (A^ + l)-electron states. Moreover, even larger clusters are needed 
to investigate the decay of the correlation effects with the distance from the local- 
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Figure 3.1: Two conformations of polyacetylene: (a) trans-PA and (b) cis-PA. 



ized extra charge. On the other hand, since rich basis sets are required to describe 
the anionic (A^ + l)-electron states properly and a powerful but rather expensive 
correlation method is employed, we reach the limits of our computational abilities 
when the size of the employed clusters reaches 30 atoms. In a bulk crystal a virtual 
WO extend in three dimensions and the number of atoms covered by the localized 
virtual orbital scales as with the size of the WO whereas in one- dimensional (ID) 
chains it only scales as Thus, for our test calculations, a linear system has been 
chosen. 



3.1.1 Polyacetylene 

A natural examples of ID periodic systems are infinitely long linear polymers. 
One of the simplest among them is polyacetylene (PA). It is a plane periodic 
molecule which consists of repeating C2H2 units and exists in two basic conforma- 
tions fFig. 13. The first one exhibits a zigzag carbon chain with alternating single 
(or long) and double (or short) C-C bonds; the C-H bonds are pointing outwards 
this chain nearly perpendicular to the molecular axis. This is trans-polyacetylene 
which is schematically shown on Fig. 13.1b . The second one is the cis conformation 
and one of its realization is drawn in Fig. 13.1b (another realization can be obtained 
by the interchange of single and double C-C bonds). 

When a polyacetylene sample is synthesized it consists of finite molecules (oligomers) 
of both conformations (at room temperature). The trans conformation is thermody- 
namically more stable and the isomerisation of the sample is achieved by heating it 
to temperatures above 150 °C for a few minutes |SI79/80|. The cis-trans conversion 
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of polyacetylene was rigorously studied in Ref. IRP+831 where the sample was an- 
nealed at 100 °C for more than 1 day. After 3 h it contained more than 50% of trans 
species and more than 80% after 1 day. The isomerization happened continuously 
during this time and homogeneously throughout the polymer rather than in isolated 
amorphous or crystalline regions. The crystallinity of these samples is estimated to 
about 90% and the amorphous content did not increase with isomerization. 

There are two methods for the synthesis of crystalline polyacetylene. In the 
work of Sirakakawa and Ikeda |SI79/80| polyacetylene was grown in the form of 
randomly oriented fibrils. Inside a particular fibril the single chains of polyacetylene 
are directed along the axis of the fibril. In order to obtain a material with parallelly 
oriented chains of polyacetylene for further crystallographic and electronic structure 
studies the grown polymer is stretched along some direction. The corresponding 
electron micrographs can be found in |SI79/80| . The characteristic diameter of fibrils 
obtained by this method ranges from 200 A to 500 A. This technique was exploited in 
a number of experimental works on polyacetylene, e.g. IFO+791 , ITG+801 , |FC+82| , 
IRP-l-831 , |YC83j and |LS-|-93| . The disadvantages of samples produced this way 
are not perfectly aligned polyacetylene chains and even the presence of a certain 
amount of randomly oriented chains, and the fact that it is impossible to investigate 
the electronic properties of tPA in the direction perpendicular to the stretching due 
to the fibrillar morphology of the material. 

The other method to produce crystalline tPA was developed at Durham Uni- 
versity f |EF80j and |EFB84j ) and provides true bulk polyacetylene. The synthesis 
utilizes a solution of a precursor polymer which is converted into polyacetylene. 
Under usual conditions, both conformations are present and the obtained poly- 
mer is amorphous. However, one can get samples with perfectly oriented chains of 
trans-polyacetylene by applying an appropriate stress at temperatures up to 120 °C. 
Then, during the conversion, an elongation of the polymer film and simultaneous 
alignment of tPA chains take place. As the result highly anisotropic free-standing 
films of 100% trans-polyacetylene are obtained. This material exhibits a compact, 
non-fibrous morphology and smooth surfaces appropriate for further optical study 
r [FL86] . |KLL87j . IBF+871 , |lJTLK87j and |L88]). 

Thus, when we refer to irans-polyacetylene as a bulk polymer we mean a crystal 
built by densely packed, flat and infinitely long tPA single chains all aligned in one 
direction. The crystalline structure of bulk tPA was studied in |FC-|-82| by x-ray 
scattering experiments and a fishbone-like packing of independent tPA chains was 
determined fFig. 13. 2j) . The space group of tPA crystals is P2i/n with a C4H4 unit 
cell consisting of two C2H2 units on neighboring tPA single chains. The parameters 
of the monoclinic structure are determined with a unique angle (3 = 91.5° and lattice 
constants a = 4.24 A, 6 = 7.32 A and c = 2.46 A. The angle between a chain plane 
and the b axis is 55°. 

As one can see in Fig. 13.21 for any atom of the tPA crystal the first few closest 
atoms belong to the same tPA chain. Therefore, the crystal consists of weakly cou- 
pled chains with strong intrachain bonding. This is evident from the high anisotropy 
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Figure 3.2: Two schematic projections of bulk trans-polyacetylene. a = 4.24 A, 
b = 7.32 A and c = 2.46 A and /3 = 91.5° (from |FC+82n . 



of the electronic properties of bulk tPA: SCF calculations give band widths between 
5 and 10 eV for crystal- momentum vectors k parallel to the chains and less than 
0.4 eV in perpendicular directions. These findings tell that the electrons in bulk 
tPA are delocalized exclusively along the polymer chain and thus each tPA single 
chain can indeed be considered as a quasi-one-dimensional crystal. Therefore, one 
can focus on the study of single ID tPA chains to describe the electronic properties 
of ^rans-polyacetylene neglecting for a moment the weak interchain interactions. 



3.1.2 Geometry of tPA single chain 

A single tPA chain consists of repeating C2H2 units with alternating short and long 
C-C bonds fFig. I3.3|) . The structure is produced by hybridization of the valence sp"^ 
electrons of carbon and the electron of hydrogen. Three out of four valence carbon 
electrons take part in this formation of a bonds (with charge density maximum 
in the molecule plane), namely one long C-C bond, one (out of two) short C-C 
bond and one C-H bond. The fourth electron occupies the 2pz atomic orbital, 
whose lobes lie out of the plane, and together with the 2pz electron of the closest 
neighbor carbon atom form the second short C-C bond which is of vr symmetry. 
The alternation in the C-C bond lengths emerges from the Peierls instability of 
ID metals |P55j . If there would be no bond alternation (sometimes also called 
" dimerization" ) and all C-C bond were equal the unit cell would consist of just 
one CH unit with one unpaired electron on the 2pz orbital per cell. In this case 
tPA would be a metal with a half-filled conduction band as can be seen from simple 
tight-binding-model considerations. However, such a ID metal can reduce its energy 
by spontaneous breaking of the translational symmetry by moving two adjacent CH 
groups towards each other and forming by this the short (or double) C-C bond. By 
this the system gains electronic energy whereas its elastic energy (i.e. the energy of 
the carbon skeleton) increases. The minimum of the system energy as a function of 
dimerization length defines the equilibrium C-C bond lengths in a (C2H2)a; chain. 
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Figure 3.3: Investigated system: irans-polyacetylene single chains. 

The reduction of translation symmetry opens the band gap and tPA becomes a 
semiconductor with a completely filled upper valence band and with a band gap 
of the order of few electron-volts. The qualitative analysis of dimerization effect in 
irans-polyacetylene is comprehensively expounded in jF95j . p. 179-183, where two 
original papers |SSH79j and [SSHSOj are summarized. 

To perform single-point ab initio calculations we need to know the positions 
of the atoms in the tPA chain. We take the relevant experimental data on the 
structure of the carbon skeleton of trans-polyacetylene from |KLL87j : the long C-C 
bond length is ri = 1.45 A, the short C-C bond length is r2 = 1.36 A and the lattice 
constant is a = 2.455 A. The angle 7 between short and long bonds is measured to 
be 122°. The measured angle does not perfectly agree with measured lengths and 
we set this angle to 121.7421° (which is well within the experimental error bar). 

The experimental data on the length and orientation of C-H bond of tPA are 
not available in the literature. Therefore, we have taken the C-H bond length 
from the paper |S92j where the geometry of tPA was optimized on both HF and 
MP2 level using different basis sets. The results for C-C bond lengths obtained on 
the MP2 level with the largest basis set (of double-zeta quality with polarization 
functions) agreed nicely with the experimental data mentioned above (ri = 1.4450 A 
and r2 = 1.3634 A) however the angle between C-C bonds was somewhat different 
(7 = 123.75°). Therefore from |S92j we have taken only the value of the C-H bond 
length which is = 1.0872 A. 

The last undefined parameter of the geometry of tPA was the angle between 
short C-C bond and C-H bond. To determine this angle we used the geometry op- 
timization option of the GAUSSIAN98 program package |C98j . We have optimized 
the angle in a hydrogen-terminated C12H14 cluster of trans-polyacetylene taking 
ri = 1.45 A, r2 = 1.36 A, rg = 1.0872 A and 7 = 121.7421° as fixed parameters. 
The angle between the two cluster-saturating C-H bonds and adjacent short C-C 
bond was set to 7 and kept fixed as well. Such an arrangement allows to reduce 
the infiuence of cluster edges on the optimized angle in the center of the cluster 
where its value is maximally close to that of the infinite chain. As follows from |S92j . 
the value of the optimized angle does almost not depend on the method (HF or 
MP2) and the basis sets employed. Thus, we performed the optimization on the HF 
level with cc-pVDZ basis sets for both hydrogen and carbon atoms |D89j . We have 
tried several starting values -^o in the range from 118° to 122° and in all cases the 
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optimized angle for the chosen fixed parameters was found to be 120.057°. 

To summarize, our test system is a plane infinitely long single chain of trans- 
polyacetylene with alternating C-C bonds fFig. 13. 3|) . There are four a bonds and 
one TT bond (on the short C-C bond) per unit cell C2H2. The lengths of long C-C, 
short C-C and C-H bonds are 1.45 A, 1.36 A and 1.087 A respectively. The angle 
between two neighbor C-C bonds is 121.742° and the angle between C-H bond and 
adjacent short C-C bond is 120.057°. 



3.1.3 Basis sets 

Throughout all our further HF and correlation calculations we use valence triple 
zeta (VTZ) atomic basis sets for carbon and hydrogen atoms |D89j (see Table ITTj) . 



These Gaussian-type-orbital (GTO) basis sets are correlation consistent, have rather 
diffuse exponents to describe properly the decay of electron wavefunctions in the 
direction perpendicular to the chain, and are rich enough to provide reasonable 
results for correlation calculations in case of the anionic (A^ + l)-electron states. 

The Gaussian-type orbitals (or more precisely Cartesian Gaussians) which were 
first proposed as an approximation to atomic orbitals in |B50j have the form 

X^{r) = Nix'y^z^e-^^"-" (3.1) 

where A'^j a normalizing prefactor and /, m and n are integer numbers which are all 
zero for s-type orbitals, one of them is equal to 1 such that they sum up to 1 for 
p-type orbitals (e.g. for 2pj., Sp^, ip^ etc. one has / = 1 and m = n = 0), they sum 
up to two (e.g. 110 or 200) for d-tjpe orbitals and so on. To approximate better 
real atomic orbitals fixed linear combinations of Gaussians are used: 



X^{r) = N^x'y^z^ i ^ c.^e-^^^^' ) (3.2) 
\j=i / 

where Cij are the so-called contraction coefficients, Qj are the exponents of the so- 
called primitive Gaussians and Li is the length of the contraction of the i-th atomic- 
like orbital. The exponents and contraction coefficients of the employed cc-pVTZ 
basis sets are listed in Table 13.11 



3.1.4 Band structure of tPA chain on the HF level 

Now, that the geometry of our test system and the basis sets are assessed, we can 
start the investigation. At first, we would like to calculate the periodic HF solution 
for an infinite tPA single chain, mean the ground-state energy (per unit cell), Bloch 
orbitals and the corresponding orbital energies (the band structure). This can be 
done by the program package CRYSTAL |(;HY98j which is able to perform HF 
calculations for periodic systems and uses Gaussian-type orbitals. However, due to 
the presence of rather diffuse Gaussians in the chosen basis set (with exponents less 
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Table 3.1: Exponents and contraction coefficients of the cc-pVTZ basis sets for 
carbon and hydrogen atoms used in all our calculations. 



atom 


type 


exponent 




contraction coefficients 










shell 1 


shell 2 


shell 3 


shell 4 


C 
















s 


8236.0000 


0.000531 


-0.000878 










1235.0000 


0.004108 


-0.000878 










280.8000 


0.021087 


-0.004540 










79.2700 


0.081853 


-0.018133 










25.5900 


0.234817 


-0.055760 










8.9970 


0.434401 


-0.126895 










3.3190 


0.346129 


-0.170352 










0.9059 


0.039378 


0.140382 


1.0 








0.3643 


-0.008983 


0.598684 










0.1285 


0.002385 


0.395389 




1.0 




P 


18.7100 




0.014031 










4.1330 




0.086866 










1.2000 




0.290216 










0.3827 




0.501008 


1.0 








0.1209 




0.343406 




1.0 




d 


1.0970 






1.0 








0.3180 








1.0 


H 
















s 


33.870000 


0.006068 












5.095000 


0.045308 












1.159000 


0.202822 












0.325800 


0.503903 


1.0 










0.102700 


0.383421 




1.0 






p 


1.407000 




1.0 










0.388000 






1.0 
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than 0.2) one faces the problem of "catastrophic" SCF results (non-converging SCF 
iterations) when default values of the ITOL parameters which control the accuracy 
of the calculation of the bielectronic Coulomb and exchange series (see |CRY98j . 
p. 74 and 123-125) are used. The effect on these parameters on the convergence of 
SCF calculations in CRYSTAL is discussed in details in |PDR88j . p. 79-85 and here 
we would only like to notice that in the case of trans-polyacetylene single chains and 
the chosen basis sets well-converged results are obtained with the ITOL parameters 
being set to 8 8 8 14 24 (the defaults are 6 6 6 6 12). The SCF ground-state energy 
per unit cell is -3.53887 eV and the iteration is converged when the energy changes 
by less than 0.14 meV (5 /iHartree) from one iteration step to the next. 

The resulting orbital energies as a function of the one-dimensional crystal-momentum 
k (the band structure) are shown in Fig. 13.41 Blue lines correspond to Bloch states 
of TT symmetry (single-electron wavefunction vanishes in the molecule plane) and 
the rest are of a symmetry. Because of the different symmetry blue and black lines 
may cross. 

From Fig. l3.4l one concludes that the ionization potential and the electron affinity 
of an infinite tPA chain are given by the top of vr-type valence band and the bottom 
of TT-type conduction band in the X point, respectively: 

IP = -£:,r(vai)(A; = vr/a), (3.3) 



EA = -£:^(cond)(A; = vr/a) (3.4) 

On the HF level these quantities are IP= 5.90 eV and EA= —0.52 eV. The band 
gap, defined as the difference of these two values, is equal to 6.42 eV. However, 
from a number of experimental studies, irans-polyacetylene is known to have a 
positive electron affinity |KP-|-89| and a band gap of approximately 2 eV |FO+79| , 
|TG-f-80| , |L88j . Therefore, trans-polyacetylene is one of many systems for which the 
HF approximation fails to give both quantitatively and qualitatively correct results. 
By inclusion of electron correlation effects, following the scheme sketched in Chapter 
2, we will approach the experimental data quite noticeably. 

In our correlation calculations we focus our attention on the five valence bands 
and the three lowest conduction bands. Taking into account general considerations 
|F95j and earlier experience of the application of wavefunction-based correlation 
methods to covalent periodic systems |BF87j . |GSF97j and |AFS00j we expect that 
the correlation effect will lead to an upward shift of the HF valence bands and a 
downward shift of the conduction bands and also to a flattening of the bands. 



3.2 Localized orbitals 

To set up the finite molecules for our correlation calculations we have chosen the fol- 
lowing hydrogen-saturated clusters of trans-polyacetylene single chain: CeHg, CgHio, 
C10H12 and C12H14. All these clusters are terminated at single C-C bonds on both 



3.2. LOCALIZED ORBITALS 



41 




Figure 3.4: The HF band structure of trans-polyacetylene. 
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edges. There dangling C-C bonds are substituted by C-H bonds having the same 
angle 7 with the adjacent double C-C bond and a characteristic C-H bond length 
of 1.087 A. Such a termination on long C-C bond exclusively preserves the flat 
geometry of clusters and also minimizes the spurital impact of the saturating hy- 
drogen atoms on the cojugated tt system of tPA. Thus, we operate with closed-shell 
molecules which perfectly reproduce some finite part of our infinite system except 
for the two terminating C-H bonds. Therefore we expect that the localized single- 
particle wavefunctions (or molecular orbitals, LMOs) in these molecules are similar 
to those in the infinite system: in particular, in the central part of the molecule 
where an LMO has the same surrounding as in the infinite system they should co- 
incides with the localized crystalline (Wannier) orbital (WO). On the edges of the 
cluster on one side there is only a terminating C-H bond which mimics the long C-C 
bond (and the remainder of the infinite chain). As LMOs centered on some bond 
always also have some small, but yet not negligible contribution of atomic orbitals 
of atoms in the neighborhood of the bond, WOs can not be reproduced correctly 
at the cluster edges. Thus, LMOs from edges of clusters will exhibit some devia- 
tions from the corresponding WOs (or equivalently from the corresponding LMOs in 
the central part of clusters). In accurate correlation calculations these edge LMOs 
are usually kept frozen or if not, the obtained results should be regarded as rough 
estimates only. 

3.2.1 The case of occupied orbitals 

To generate localized occupied spatial SCF orbitals (i.e. the bonding orbitals) on 
the clusters we use the Foster-Boys localization scheme |FB60j which is implemented 
in MOLPRO |M2000j . The idea of this scheme is to find that unitary matrix Uij 
in Eq. (j2.26|) which minimizes the sum of the spread of Norh pre-selected spatial 
orbitals: 

1=1 

where the orbital spread is defined as 

crl = y"(r-R.)^|(/..(r)|W, 
and the center of z-th localized orbital ipi as 

K, = j r\^,{r)\'dV (3.7) 

The question arises which canonical orbitals we have to provide for the localiza- 
tion to get properly localized occupied LMOs, which correspond to the individual 
bonds. We have count the valence electrons in the molecule (i.e. the electrons on 
the outer shells of each atoms which take part in the formation of bonds) and half 
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Figure 3.5: Three different a bonds in irans-polyacetylene clusters: (a) C-H bond, 
(b) long C-C bond and (c) short C-C bond. The contour values increase in steps 
of 0.025 a.u. 




Figure 3.6: One of the vr bonds of a CgHio cluster of tPA. Contours are in steps of 
0.025 a.u. The plot plane is parallel to the molecule plane at a hight oi z = 0.5 A. 

of this number defines the number of the occupied canonical orbitals to be local- 
ized starting with the HOMO (highest occupied molecular orbital) and continuing 
with the successive energetically highest orbitals. All the other occupied orbitals 
are core orbitals which do not participate in bonding. In our case they are linear 
combinations of Is orbitals of carbon atoms. 

A special remark has to be made for flat molecules where orbitals of two different 
symmetries {a and vr) exist. Canonical orbitals of vr type can be readily found in 
the set of SCF orbitals of a flat molecule as they have zero value in the molecule 
plane. Then the canonical a and vr orbitals can be localized separately and LMOs 
of a and ir symmetries are obtained. By doing so one reduces the number of LMEs 
to be treated on the correlation level as off-diagonal LMEs between a and tt orbitals 
are zeros and the effort for the correlation calculations, which are the most difficult 
and time consuming part of the whole approach, is reduced substentially. 

In Fig. 13.51 contour plots of localized one-particle wavefunctions </3a(r) for three 
different a bonds of a CgHio cluster are shown: (a) C-H bond, (b) long C-C bond 
and (c) short C-C bond. Red color denotes positive values of V5a(i") and blue color 
denotes negative ones. 
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Figure 3.7: C10H12 cluster of tPA with the names given to each bond. The reference 
(central) unit cell is indicated by the dashed frame. 

In Fig. Klfil a typical localized vr orbital is shown. As it vanishes in the molecule 
plane {ip{r)^z=o = 0) the plot plane for the contour plot is shifted to z = 0.5 A. 
All plots in Figs. 13.51 and 13.61 are obtained by the MOLDEN program package 
IMOLDENj . 

Comparing Fig. 13.51 and 13.61 we see that the vr-bonds are more diffuse than 
the (T-bonds. The former have also a quite noticeable contribution of atomic basis 
functions from carbon atoms in the two nearest-neighbor unit cells. Therefore, for 
an accurate treatment of the correlation corrections to LMEs corresponding to vr 
bonds substantially larger clusters are needed than for the case of a bonds. 

Having obtained the occupied LMOs we calculate the local SCF matrix elements 
IPab as defined in ()2.46j) and identify them with the corresponding crystalline local 
matrix elements IPn^nn' which enter the truncated sum ()2.32|) . The hopping matrix 
elements decay with the distance between two WOs. Setting some threshold for the 
absolute values of LMEs in ()2.32|) we define the two most distant WOs and by this 
the minimal size of the cluster in which the LMEs are calculated. Also, for choosing 
the proper cluster one takes into account that the LMOs on the edges of clusters 
can not reproduce the corresponding WOs correctly. Thus, the size of the cluster 
on which the LMEs are calculated is given by the two most distant WOs which give 
hopping matrix elements above the threshold plus one additional unit cell on both 
edges of the cluster. In the case of tPA the threshold is set to 1 mHartree (0.027 
eV) , and the minimal-sized saturated cluster arising from that is a C10H12 molecule. 

This cluster is schematically shown in Fig. 13.71 For convenience we give names 
to the different bonds in the unit cell: on the double C-C bond we distinguish bonds 
of TT and a type and call them "pi" and "sig", respectively, the single C-C bond is 
called "CC" and the two C-H bonds are labelled "CHI" and "CH2". In the central 
unit cell which is marked by dashed rectangle the bond names are taken as they 
are. Bonds in the other unit cells acquire an additional index which correspond to 
the ordering numbers of the cell they are located in starting from the central 
(reference) cell and being negative on the left- and positive on the right-hand side. 
We build the finite matrix IP_R,„n' between the WOs ipon and ipun' defined in ()2.30p 
with R = Ucai. The SCF values of IP R,nn' matrix are summarized in Table 13.21 
for (T-type WOs and in Table 13.31 for vr-type WOs. The number of relevant non- 
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Table 3.2: SCF values of the IPi?,,nn' matrix for cr-type WOs (in eV). 





R 


-n.,CHi n' 




1 1^ T~) /^TTri / 

-ft,CH2 n' 


IP-R,CC n' 
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-0 696 
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-0 767 


CHI 


n 


18 Qn4 


3.022 


-0 856 


81 7 




n 


3.022 


22 396 


3.022 


3.022 


CH2 





-0.856 


3.022 


18.904 


2.968 


CC 





0.817 


3.022 


2.968 


21.230 


CHI 


-lai 


0.338 


0.634 


-0.803 


2.968 


sig 


-lai 


-0.190 


-0.696 


0.634 


3.022 


CH2 


-lai 




-0.190 


0.338 


0.817 


CC 


-lai 




0.202 


-0.153 


-0.767 


sig 


-2ai 




-0.049 




0.202 


CHI 


-2ai 








-0.153 



Table 3.3: SCF values of the IP R^nn' matrix for vr-type WOs (in eV). 



IPo.pi pi 


IP±a 


pi pi 


TP TP TP 

±2a,pi pi ±3a,pi pi ±4a,pi pi 


10.526 


1.6 


B3 


-0.365 0.120 -0.047 



equivalent entries in this matrix is finite since we omit hopping matrix elements 
which are below the threshold. The matrix IP/j „„/ is clearly diagonal dominant 
and the matrix elements nicely decay with increasing distance of the outers of the 
respective Wannier orbitals. 

Let us emphasize that the matrix IP n^nn' contains orbital energies and hop- 
ping matrix elements for an infinite tPA chain while its entries are calculated on the 
C10H12 molecule being IPafe. To build the matrix IP R^nn' correctly we use both, trans- 
lation and inversion-center symmetry, of the infinite system and usually those matrix 
elements of IPa^ which correspond to bonds close to the center of the molecule. For 
example, all diagonal matrix elements IPonn are associated with IPaa from the central 
unit cell. For off-diagonal elements we have IP2a^si^ sig = IP^2Xsig sig = IP™ri,-i), (sig,i) 
(use of translation symmetry) and IPo|'ccch2 = IP-^a,cc chi = IPcc, ch2 (use of 
inversion-center symmetry). Here, for the distinctness we temporarily use super- 
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Figure 3.8: Valence bands of irans-polyacetylene obtained on the SCF level by LMEs 
calculated on clusters (red) and by a CRYSTAL calculation on an infinite tPA single 
chain (black). 



scripts "chain" for LMEs from the infinite chain and "mol" for LMEs obtained from 
SCF calculations of molecules. 

Substituting the matrix elements from Tables 13.21 and 13.31 in the truncated sum 
(|2.32p one obtains Nj. matrices IP„,„/(/Cj) on the chosen /c-mesh of the Brillouin zone 
(z = 1, 2, . . . , Nk). After diagonalization these matrices provide the positions of 
the bands at each point ki. By choosing a rather dense grid one can easily produce 
a smooth band structure. In Fig. 13.81 the five valence bands of trans-polyacetylene 
obtained this way are presented and compared with the SCF bands of an infinite 
tPA single chain already shown in Fig. 13.41 Comparison reveals that the LMEs 
reproduce the SCF band structure very nicely and only small deviations (less than 
0.2 eV) from the band structure by CRYSTAL appear due to the truncation in 
Eq. ()2.32j) . This demonstrates that LMEs obtained in molecules are practically 
equal to those comming directly for the infinite system using Wannier orbitals. 
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3.2.2 The case of virtual orbitals 

The situation becomes more difficult when locahzed virtual orbitals have to be gen- 
erated. In general, one can not simply use the same strategy as for the localization 
of the occupied orbitals. Occupied LMOs are obtained by applying e.g. the Foster- 
Boys localization scheme to the A^bond energetically highest occupied canonical spa- 
tial orbitals where iVbond is the number of bonds in the molecule. If one applies the 
localization procedure to the A^bond energetically lowest unoccupied canonical spatial 
orbitals one usually does not get any localized virtual orbitals (or antibonds) which 
can be regarded as the unoccupied WOs corresponding to the lowest few conduction 
bands. This difference in getting bonds but no antibonds lies in the fact that the va- 
lence bands are well-separated from the rest of the bandstructure by the band gap 
and that accordingly the occupied spectrum of the molecule representing a finite 
part of the crystal is clearly separated from the virtual spectrum by the HOMO- 
LUMO gap (see Fig. 13. 9|) . In contrast, the lowest conduction bands are usually 
entangled with higher conduction bands and therefore the spectra of virtual states 
corresponding to lower and higher conduction bands overlap. For instance, in the 
band structure of tPA the five lowest conduction bands (which would correspond to 
the five antibonds per unit cell) are not separated from the other conduction bands. 
In this case one can not select in a simple and controlled way some suitable A'^bond 
unoccupied canonical orbitals which could be converted to antibonds by the local- 
ization procedure. To illustrate this we show on Fig. 13.91 the spectrum of a CgHg 
cluster together with the SCF band structure of tPA. 

Alternatively, one could localize all virtual canonical orbitals of the molecule. 
However, even the energetically lowest localized virtual orbitals obtained this way 
have too high orbital energies since they contain contributions of both energetically 
low- and high-lying canonical orbitals. In this case non of the virtual LMOs can 
be regarded as a WO corresponding to the lowest few conduction bands. Only the 
whole set of LMOs can be used to reproduce all conduction bands. Therefore LMOs 
obtained by localization of all virtual molecular orbitals can not be used in our 
approach for correlation correction to the few lowest conduction bands. 

Thus, in general, antibonds can not be obtained on clusters by localization of 
some set of virtual canonical molecular orbitals (except of the cases of well separated 
conduction bands). Hence, we propose a completely different way of getting localized 
orbitals in clusters which can reproduce the bands of our interest: Wannier orbitals 
obtained in an infinite system are projected onto a cluster. If the Wannier orbital 
can be obtained as a linear combination of atomic orbitals with contributions of 
the latter rapidly decaying with the distance from the orbital center (see Fig. I3.1(jj) 
such a WO can be represented on the cluster when the contributions from those 
atoms which are not present in the cluster are negligibly small such that they can 
be omitted. Then the projected WO can be used as the LMO. For instance, the 
WO shown schematically in Fig. 13.101 can reasonably be represented in a cluster 
consisting of five unit cells or more. 
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Figure 3.9: SCF band structure of tPA (on the left) and the one-particle spectrum 
of a CeHg cluster (on the right). 
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-3a -2a -a a 2a 3a 

Figure 3.10: Schematic diagram showing the contributions of the atomic orbitals 
from different unit cells to a WO centered on the zero unit cell of a ID crystal. 
Horizontal bars denote the unit cells and the hight of the vertical bars is proportional 
to the squared norm of the wavefunction contribution (or population) from the 
individual cells R = 0,±a, ±2a, .... 



So, we have first to get well-localized Wannier orbitals in the infinite system 
which correspond to the conduction bands of interest. The concept of maximally lo- 
calized Wannier orbitals was introduced for an isolated band or composite bands (an 
isolated group of bands) by Marzari and Vanderbilt |MV97j . There, the molecular 
localization criteria of Foster and Boys were applied to minimize the total spread of 
the WOs in the periodic system. In |MV97j Bloch functions ^/'nk(r) represented on a 
discrete mesh of /c-points (n labels the band) are used and the Wannier orbital in cell 
R associated with band n is determined as w„(r — R) = 1/Nk e~*''''^'?/'nk(r) {Nk 
is the number of considered k points in the Brillouin zone). The minimization of the 
total spread functional is carried out in the space of unitary matrices Uml describing 
the rotation among the Bloch functions at each fc-point: ip'^^^ = Umn'ipm]i- 

Later this idea was generalized in |SMV01j for the case of entangled bands to ac- 
count for the case where no separate group of bands can be found. There, one speci- 
fies an energy window with A^tot(k) bands with the Nb bands on interest being present 
together with higher or/and lower lying bands (at each /c-point A'tot(k) > Nh). At 
fc-points where A'tot(k) > Nb an iVfe- dimensional subspace of Bloch states which 
correspond to the bands of interest is first established (by minimization of some 
functional) and then the localization method from |MV97j is applied to the selected 
Nb states at each A;-point. This way the maximally localized Wannier orbitals asso- 
ciated with selected bands are obtained. 

A variant of the method of |MV97j for obtaining well-localized Wannier orbitals 
corresponding to composite bands was recently implemented in the CRYSTAL code 
( IBZ+Ol] and jZWDSOT] ). The Bloch orbital for band n at point k is represented in 



CRYSTAL as the linear combination of atomic orbitals x^(r — s^) (being contracted 
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Gaussian-type function) centered at an atomic site 

M 

^nk(r) = ^a;:(k) 5^e^^-^x.(r - s, - R) (3.8) 

ti=l R 

where n labels atomic orbitals centered in the reference cell and R denotes the lattice 
vector. The Nf, Wannier orbitals associated with the A'^f, composite bands at each 
unit cell are then obtained as linear combinations of atomic orbitals 

M 

Wr.{r-K) = Y,Y1 - - (3-9) 

^l=l R' 

with the orthonormality condition being fulfilled {wn{r — lV)\wn' (r — R')) = 5nn' 5r,r' 
This new feature of CRYSTAL was explicitely used in our work to get well-localized 
virtual WOs in trans-polyacetylene corresponding to the three lowest conduction 
bands since they can be regarded as composite bands being separated from other 
conduction bands by a gap in the whole Brillouin zone (see Fig. l3.4l or l3.9p . However, 
at this point we would like to note that the idea of band disentanglement originally 
developed in |SMV01j was very recently implemented in the CRYSTAL program 
package by our work group and WOs corresponding to the lowest conduction bands 
of bulk diamond and silicon were obtained |IBj . Thus, the limitation of the method 
reported in |BZ+01| and [TWDSOlj for separated group of bands is removed and 
WOs for any bands of interest can now be obtained for the general case. 

Having obtained well-localized virtual WOs in periodic system as linear combi- 
nations of atomic orbitals centered on atoms of a finite region of the crystal (usually 
some few unit cells close to the center of the WO) we need to project these orbitals 
onto the basis set of a finite molecule being a hydrogen-saturated cluster of the 
crystal. Then we could use them as virtual LMOs to get LMEs in our approach and 
also for further correlation calculations. For this purpose the program CRYREAD 
CRYREADj was recently designed in our group. It reads the WOs obtained in 
CRYSTAL as linear combinations of atomic functions and projects them onto the 
basis set of the user-specified cluster of the periodic system. The program takes all 
those WOs with norms exceeding a user-specified threshold after the coefficients c^^: 
at atomic orbitals Xfi{^ ~ ~ -^') centered outside the cluster are set to zero, makes 
them othonormal and writes the corresponding matrix of coefficients c^'^/ where 
vectors R' and R are restricted to the specified cluster. After some further mod- 
ification (contraction coefficients of atomic orbitals centered on cluster-saturating 
hydrogen atoms are added and set to zero) this matrix can be read by the program 
package MOLPRO and virtual localized molecular orbitals are obtained directly in 
the form of linear combinations of basis functions of the molecule. Virtual LMOs 
obtained this way perfectly represent WOs of the periodic system. 

To be more precise, in the case of irans-polyacetylene we generate with CRYS- 
TAL Wannier orbitals associated with core, valence and the first three conduction 
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bands, performing three separate calculations for each subset of bands. The pro- 
cedure for getting well- localized crystalline orbitals in CRYSTAL ( IBZ+Ol] and 
|ZWDSOT] ) is an iterative one and requires (especially for virtual orbitals) an initial 
guess for the Wannier orbitals to start with. As discussed in |ZWDSnT] the choice 
of the initial guess is critical for systems with fully covalent bonds (like diamond, 
silicon and also irans-polyacetylene). Therefore, we have to design good starting 
Wannier orbitals for each subset of bands. For core bands Is atomic orbitals of car- 
bon atoms are provided. In the case of valence bands we expect that the associated 
WOs are similar to the molecular bonds of hydrogen-saturated clusters. Therefore, 
we specify five bonding pairs of atomic orbitals per unit cell which correspond to 
the five bonds of the covalent system. The choice of the initial guess for three un- 
occupied WOs associated with the three lowest conduction bands of tPA is more 
subtle. We suppose that these three WOs correspond to three out of five antibonds 
per unit cell. As one of the conduction bands is formed by states of vr symmetry the 
corresponding WO should be the vr-antibond centered on the short C-C bond. For 
starting WOs associated with the two cr-bands the C-H antibonds are used. This 
choice is based on the fact that the C-H bonds in hydrogen-saturated clusters of 
tPA are energetically higher than the C-C a bonds and therefore the C-H antibonds 
are expected to be energetically lower than the C-C a antibonds. Since we already 
have occupied WOs Wn{r) (or bonds) the antibonds of interest are obtained by a 
change of the sign of the coefficients c^'^ from the atomic orbitals Xaj(i" — s^) centered 
at one of two atoms of the corresponding bond. Providing antibonds constructed 
this way we were able to generate well-localized unoccupied WOs associated with 
three lowest conduction bands of tPA. The condition of orthonormality is fulfilled 
for all Wannier orbitals: (wi(r — R)|wj(r — R')) = ^jj^R^R' (indices i and j running 
over core, valence and the first three conduction bands). 

By the program CRYREAD we project the virtual WOs onto a C12H12 clus- 
ter (six unit cells, no saturating hydrogens) of the infinite tPA chain. Thus we 
obtain 18 (3x6) WOs centered on this cluster. The coefficient matrix produced 
by CRYREAD is enlarged to be consistent with the set of basis functions of an 
isolated C12H14 molecule (hydrogen-saturated). The enlarged matrix is read by the 
MOLPRO program where the SCF calculation on the C12H14 molecule has been per- 
formed and the set of one-particle SCF wavefunctions have been obtained. Then, 
the occupied and virtual spaces in C12H14 are constructed as follows. The occupied 
space {ipa} is taken from the SCF calculation on C12H14 and is not changed. Next, 
the WOs Wnir — R) = <^r(r) (obtained in the periodic system and projected onto 
the C12H12 cluster) are provided and are made orthogonal to the occupied space of 
the C12H14 molecule by Schmidt orthogonalization and re-orthonormalized there- 
after by Lowdin (or symmetric) orthogonalization (see |S096j . p. 142-145). After 
this orthogonalization step the obtained virtual LMOs ip'^ may slightly differ from 
the projected WOs ipr as the occupied space of the C12H14 molecule differs from the 
one of the embedded C12H12 cluster. The space of the remaining virtual one-particle 
orbitals of C12H14 {ip'r} is constructed to be orthogonal to the occupied space {^a} 
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(b) 




Figure 3.11: One C-H (a) and one C-C tt (b) antibonds of a C12H14 cluster of tPA: 
the plane for the contour plot of the C-H antibond lies in the molecular plane, that 
of the C-C vr antibond is lifted upwards to z = 0.5 A. Red lines correspond to 
positive and blue lines to negative vlues of V5r(r)- The contour values increase in 
steps of 0.025 a.u. 



and the specified subspace of the 18 virtual LMOs {v^^}- 

The 18 virtual LMOs constructed this way turned out to be 6 C-C antibonds of 
TT type and 12 C-H antibonds of the C12H14 molecule. Two of these antibonds (one 
C-H and one C-C) from the central part of the molecule are presented in Fig. 13.111 
Again as vr orbitals have zero values in the molecule plane the plot plane in the 
contour plot in Fig. I3.1ir b) is lifted by 0.5 A. As one sees from these plots, the 
unoccupied LMOs are substantially more diffuse than the corresponding occupied 
ones fFig. 13.51 and 13. 6|) . 

So, we got virtual LMOs of the molecule ip'^. in form of linear combinations of basis 
functions Xfj. centered on atoms of this molecule. The LMOs represent WOs and can 
be used to obtain LMEs on the SCF level and to reproduce the three conduction 
bands. However, at this point we would like to turn back to the question of getting 
virtual LMOs by applying the Foster-Boys localization scheme to virtual canonical 
orbitals ipr- The problem was that, in general, one can not simply select a subset of 
consecutive virtual canonical orbitals. Now with the LMOs at hand we can project 
the virtual LMOs onto the space of the virtual canonical orbitals {ipr} of the C12H14 



3.2. LOCALIZED ORBITALS 



53 



molecule (which altogether consists of 359 orbitals): 

359 

\v's) = J2\'^r){^PM) (3.10) 

r=l 

and select the required subset for the localization by chosing those with the largest 
projection coefficients {iprl'^'^)- This procedure is justified as long as all other co- 
efficients are close to zero which was indeed the case for all our clusters. In case 
of the hydrogen-saturated C12H14 cluster of irans-polyacetylene the 20 energetically 
lowest virtual canonical orbitals were selected by the projection ()3.10p . After local- 
ization they give 6 C-C antibonds of vr symmetry and 14 C-H antibonds including 
the two C-H antibonds with saturating hydrogen atoms which are not present in 
the crystal. (Again canonical orbitals of vr and a symmetry are localized separately 
from each other.) The LMOs obtained this way will be denoted by ipr- We have 
to emphasize here that we could get these LMOs ipr by Foster-Boys localization of 
suitably selected virtual canonical orbitals only because of two unique properties of 
irans-polyacetylene: i) it has three separated conduction bands (see Fig. 13. 4|) which 
correspond to the C-C vr antibond and two C-H antibonds (note that no separated 
bands corresponding to two C-C a antibonds can be found in the band structure 
of tPA) and ii) the two C-H bonds which connect two carbon atoms on the edges 
of the cluster with the saturating hydrogen atoms are similar to the characteristic 
C-H antibonds of the (C2H2)2: ID crystal and the energy of an electron in the ar- 
tificial antibond is approximately the same as the one in the C-H antibond of the 
crystal. In general, localization of virtual canonical orbitals fails to produce proper 
antibonds. For instance, no suitable subset can be found in the virtual spectrum 
of hydrogen-saturated clusters of diamond (or silicon) to get the C-C (or Si-Si) 
antibonds by Foster-Boys localization. Thus, in general, one has to use projected 
WOs as LMOs. We decide to use the LMOs ipr for irans-polyacetylene because they 
are easier to handle than the projected WOs ip'^ but we would like to analyze first 
the difference between and (p'^.. 

For the C-C vr antibonds there is no noticeable difference between the projected 
WOs and the LMOs obtained by localization of the first 6 virtual canonical vr orbitals 
as there is no effect of the saturating C-H antibonds of a symmetry on the C-C vr 
antibonds. Therefore, the LMEs corresponding to states with an attached electron 
in a projected virtual WOs and an LMOs of vr symmetry are practically the same. 
In Table Owe list the LMEs EAR^rnm' defined in (imH) for {N + l)-electron states 
where the additional electron is placed in a vr antibonds as obtained by Foster- 
Boys localization. The names of antibonds coincide with those of the bonds but 
are marked with an asterisk. Again a rapid decay of the matrics elements with 
increasing bond distance is clearly discernable. 

When C-H antibonds are obtained by localization of the 14 energetically lowest 
unoccupied canonical cr orbitals of C12H14 with the latter ones containing contibu- 
tions of the two C-H antibonds with the saturating hydrogen atoms. Since unoc- 
cupied LMOs are more diffuse than occupied ones (and spread over three unit cells 
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Table 3.4: SCF values of the EA^ ^m' matrix for vr-type unoccupied WOs as obtained 
on a C12H14 cluster (in eV). 



EAo,pi* pi* EA±a^pi* pi* EA±2a,pi* pi* EA-|-3a^pi* pi* EA-|-4a^pi* pi* 

-4.497 -1.054 0.701 -0.146 0.025 
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Figure 3.12: Conduction bands of trans-polyacetylene obtained on the SCF level by 
LMEs calculated in a C12H14 molecule with projected unoccupied Wannier orbitals 
(red) and by CRYSTAL calculation on the infinite tPA single chain (black). 

of tPA as one sees in Fig. I3.11|) the C-H antibonds obtained in the molecule are 
affected by the presence of the saturating atoms. As a consequence, the local SCF 
matrix elements for attached electron states EA^ mm' of a type calculated in the 
C12H14 molecule by means of the projected WOs as localized virtual orbitals differ 
slightly from those calculated with the localized orbitals in the molecule. The two 
relevant matrices are listed in Table 13.51 The first two colunms of data show the 
values of the LMEs obtained when projected virtual WOs of a type are used as C-H 
antibonds. The last two columns represent the same matrix elements obtained with 
localized virtual orbitals from a Foster-Boys localization of the canonical a orbitals 
of C12H14. They differ by up to 0.12 eV which is small but not negligible for the 
calculation of the SCF bandstructure. The difference between the matrix elements 
EA/{ is caused by the presense of the two saturating hydrogen atoms in C12H14 
molecule. 

To demostrate how the corresponding virtual LMEs reproduce the conduction 
bands we substitute the data from Table 13.41 and the first two columns of Table 13.51 
into the truncated sum (j2.33|) and compare in Fig. l3.12l the bands obtained this way 
(red lines) with the three first conduction bands of infinite tPA chains as obtained by 
CRYSTAL (black lines). The deviation of bands by LMEs from those by CRYSTAL 
does not exceed 0.2 eV for all /c-points like in the case of the valence bands and is 
mainly due to the truncation of the infinite summation (j2.33|) . 

It has to be noticed that one needs hopping matrix elements between more dis- 
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Table 3.5: SCF values of the EA^ ^m' matrix for C-H antibonds in C12H14 cluster 
(in eV). 







Projected WOs 


LMOs by FB localization 


m' 


R 


EA/j^CHl* m' 




EAr^CHI* m' 


EAr^CH2* m' 


CH2* 


Aai 


-0.063 




-0.068 




CHI* 


3ai 


0.120 


-0.083 


0.092 


-0.054 


CH2* 


3ai 


0.097 


0.120 


0.107 


0.092 


CHI* 


2ai 


-0.289 


0.126 


-0.244 


0.078 


CH2* 


2ai 


-0.189 


-0.289 


-0.177 


-0.244 


CHI* 


lai 


1.012 


-0.236 


0.914 


-0.194 


CH2* 


lai 


-0.031 


1.012 


-0.025 


0.914 


CHI* 





-5.603 


-0.013 


-5.487 


-0.042 


CH2* 





-0.013 


-5.603 


-0.042 


-5.487 


CHI* 


-la; 


1.012 


-0.031 


0.914 


-0.025 


CH2* 


-la; 


-0.236 


1.012 


-0.194 


0.914 


cm* 


-2ai 


-0.289 


-0.189 


-0.244 


-0.177 


CH2* 


-2ai 


0.126 


-0.289 


0.078 


-0.244 


CHI* 


-3ai 


0.120 


0.097 


0.092 


0.107 


CH2* 


-3ai 


-0.083 


0.120 


-0.054 


0.092 


cm* 






-0.063 




-0.068 
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Figure 3.13: Two lowest conduction bands of a-symmetry obtained on the SCF level 
by LMEs calculated fully in C12H14 molecule (red lines) as compared to a bands 
from CRYSTAL calculation on the infinite tPA single chain (black lines). 

tant virtual LMOs (up to 4-th nearest neighbor cell) than in the case of the occupied 
LMOs to reproduce the a conduction bands with the same accuracy as valence bands 
(where 2nd nearest neighbor cells were enough). Comparing Tables and IT^ we 
see that the decay of the SCF hopping matrix elements with the distance between 
the localized virtual orbitals is slower than in the case of the occupied orbitals since 
the former are more diffuse (cf. Fig. I3.5f a) and l3.HT a)). This implies that one 
needs larger clusters for getting the EA^^mm' matrix elements than the IPn^nn' ma- 
trix elements if the same value for the truncation threshold in (j2.33|) and ()2.32|) is 
used. 

In the case of irans-polyacetylene we have the unique possibility to compare the 
local matrix elements EA/j ^m' obtained with projected unoccupied Wannier orbitals 
(f'^ and antibonds (pr obtained on a cluster by Foster-Boys localization of canonical 
virtual orbitals ^/^^ of this cluster. If we also produce the lowest a conduction bands 
by substitution the matrix elements from the two last columns of Table 13.51 into 
fl2.33|) and compare them to the corresponding a bands obtained by CRYSTAL for 
the infinite tPA chain (see Fig. I3.13j) we see rather big deviation for /c-points close 
to the X point. This is most likely due to the onset of an avoided crossing in the 
upper a conduction band close to the X point (see Fig. I3.4p . Therefore, for the SCF 
values of the EA^^mm' matrix elements we take those which are obtained with the 
projected unoccupied WOs. 

3.3 Correlation corrections to the local matrix el- 
ements 

The aim of our approach is to obtain the local matrix elements as defined in (j2.63p 
and ()2.64|) with the correlation effects taken into account. This is equivalent to 
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finding corrections 6H^°" and 6Eq°" to the matrix elements Hij from Eq. ()2.74j) and 
tlie ground-state energy Eq obtained in tlie framework of the HF approximation, 
respectively: 

SH^°" = H^" - H,, = H^- Hij (3.11) 

and 

SE^°" = El"'' -Eq. (3.12) 

Then, the local matrix elements including the correction effects are obtained by 
adding these corrections to the corresponding SCF matrix elements as follows from 
Eqs. (ESni), (ESI, (nnT|l and (1^1^ : 

IPr = IPafe + ^IPar = Hat + SHT - SabiEo + ^E^) (3.13) 

and 



EA^r = EA,, + 5EA^r = -{Hrs + SH^:") + 5,,(Eo + 6El"") (3.14) 
The correlation corrections to the local SCF matrix elements 

61FT = - Sab6Er (3.15) 

and 

5EA^°"- = -6H',^'' + 6rs6El°" (3.16) 

appear as the consequence of electron correlation effects in the charged and the 
neutral system (which are two different effects as explained in chapter 2.4). Formally, 
SW^f^' and 6EA'^°J' are non-zero because of the CI expansion of the wavefunctions 
according to Eqs. (j2.76p and (j2.68|) . Without this expansion (i.e. when SCF many- 
particle wavefunctions are used) we get the SCF local matrix elements IPab and 
EArs which are listed in Tables IT^3. 51 

In principle, correlation corrections to all matrix elements listed in these tables 
can be calculated. However, in order to get correlated band structure with sufficiant 
accuracy one needs to account for correlation effects only on the diagonal and the 
few largest off-diagonal matrix elements. As known from earlier works |GSF93j . 
GSF97j . |ARM98j and |AFS00j (and as confirmed by the present work) the corre- 



lation corrections SIP^^" and 5EA™" are small fractions of the corresponding SCF 
matrix elements IP at and EA^^. Also these corrections are obtained with final ac- 
curacy (10% say). Therefore one can conclude that for a given set of SCF local 
matrix elements reproducing some set of bands on the SCF level there is no sense to 
calculate correlation corrections to LMEs which are of the order of (or smaller than) 
10% of the largest matrix element in the given set since the obtained correction to 
such small LMEs i) has no noticeable effect on the band structure and ii) themselves 
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are of the order of the error bar for the largest matrix element. Thus we can restrict 
the number of matrix elements which need to be corrected to all diagonal elements 
and all those off-diagonal elements whose absolute values exceeds 10% of the largest 
diagonal matrix element in the given set (see Tables IH.2HH.5|1 . 

Practically, these off-diagonal elements are those where the two involved localized 
orbitals have the largest overlap. So, for n orbitals (occupied and virtual) the 
hopping matrix element for two vr bonds (or antibonds) from two adjasent unit 
cells are considered. In the case of the occupied a bonds we calculate correlation 
corrections to off-diagonal elements ^IPa™"^ for all pairs of bonds sharing one atom. In 
tPA there are three such non-equivalent hopping matrix elements between adjacent a 
bonds: IPcH2,sig5 IPch2,cc and IPcc.sig (see Fig. l3.7p . The C-H antibond is extended 
along the chain and therefore has the largest overlap with its translational copy 
from the nearest-neighbor unit cell (see Fig. IH.llk ). Thus, taking into account the 
inversion and translation symmetry of tPA we only need to care for one off-diagonal 
matrix element, EAch2*,ch2,i*, to correct the a conduction bands. 

As discussed in chapter 2.4, electron correlation is mainly a local effect. It takes 
into account the possibility of electrons in two neighbor localized orbitals to avoid 
each other and, by this, to reduce the energy of electron-electron repulsion. Also the 
rearrangement of electrons due to the presence of an extra charge happens mainly in 
the close vicinity of the charge. Therefore, the major contribution to the corrections 
^IP^^"^"^ and 6EA'^°" to the local SCF matrix elements IPa^ and EA^s is produced by 
correlation of electrons in localized orbitals which are located in unit cells close to 
orbitals a and b (or r and s respectively). Thus, for each matrix element IPab (or 
EArs) one can explicitly define a cluster in which electron correlation gives the major 
contribution to ^IP^""^ (^EA™"). For this cluster correlation should be treated with 
high accuracy (that is why the multireference CI method is used in our work). 

To account for correlation effect on IPab matrix elements coming from electrons in 
(J bonds a CgHg cluster is enough. However, in order to treat accurately correlation 
between electrons in tt bonds one needs CgHio and C10H12 molecules because the 
occupied n bonds are more diffuse than the a bonds. As unoccupied localized 
orbitals are even more diffuse than the occupied ones, correlation calculations in 
C12H14 molecule have to be performed to properly account for corrections to the 
EArs matrix elements. The details of the application of the multireference single and 
double configuration interaction method to get the corrections ^IP^b" and ^EA™" 
are outlined in the next section. 

3.4 The method of local increments 

Having chosen the molecule (a hydrogen-saturated cluster of the crystal) for explicit 
correlation calculations we have to treat all electrons of the (A^ — 1), or (A^ + 1)- 
electron system on the correlation level except of the electrons in the core orbitals 
and in the bonds with the saturating hydrogen atoms. Let us say that in some cluster 
we have to handle electrons from nocc occupied bonds. Then, to get the correlation 
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correction to the ground-state energy all singly and doubly excited configurations 
with electrons removed from these riocc bonds are included into the expansion ()2.68|) . 
For hole states, Uqcc reference configurations \Ef?°'"'" are formed as linear combinations 
of TT-occ one-particle configurations $a with the hole in one of the tt-occ bonds. Then, 
excited configurations are produced by taking one or two out of (2nocc — 1) electrons 
from the occupied bonds and placing them into some virtual orbitals as defined in 
Eq. ()2.76p . and the energies of nocc hole states have to be calculated. In the case of 
an attached electron in n^ij-t antibonds of the cluster we have nvirt reference states 
with the extra electron being delocalized over these antibonds. Due to the presence 
of the extra electron the electrons in the Uqcc occupied bonds rearrange and in order 
to account for the correlation effect in this case all excited configurations with one or 
two electrons removed from either one of the nocc occupied or virtual orbital which 
carries the extra electron and placed in some free virtual orbitals have to be taken 
into account in Eq. ()2.7(jj) . Then, the energies of nvirt states are calculated. In all 
these three cases the number of excited configurations which are taken into account 
is huge even for small molecules and robust CI calculations as described above are 
not feasible (an enormous number of coefficients a from Eqs. ()2.68|) or ()2.76|) have to 
be optimized simultaneously). To proceed one needs an approximate scheme which 
consists of a series of calculations with the number of excited configurations being 
drastically reduced but which still accounts for electron correlation in the cluster 
with the desired accuracy. 

3.4.1 Formulation in terms of bonds 

A simple but powerful scheme to estimate correlation correction to SCF ground- 
state energy of clusters and solids was proposed by StoU in Ref. !St92aj and |St92bj . 
This scheme is based on a clear physical idea and systematically sums contributions 
to 5Eq°^^ where each term accounts for correlation of a very limited number of 
electrons in localized orbitals. The method is as follows. First, corrections to the 
ground-state energy from electrons in each bond a individually are calculated: = 
i?Q°''''(a) — £"0, the energy E'^"{a) is obtained from Eq. (|2.5ip where in the CI 
expansion for the correlated wavefunction (j2.68|) only electrons in bond a are used 
to construct the excited configurations and Here a stands for a spin orbital 
with the same spatial orbital as a but with opposite spin-orientation with respect 
to the spin orbital a. Let us notice here that a drastic reduction of the number of 
excited configurations in (j2.68|) takes place when excitations of two electrons only 
are considered. Summing up the contributions from all bonds one gets a first 
estimate 5E'^^ for the correction to the HF ground-state energy. To account for the 
non-additive corrections ea,b = EQ°"{a,b) — (e^ + e?,) — Eq, correlation calculations 
with only the pair (a, b) of bonds being open for producing excited configurations 
are performed. Adding all pair corrections ea^b to SE^^^ we reach the next order of 
approximation 6Eq . The third order is obtained when we account for correlation 
corrections of electrons in bond triples (a, 6, c) in the cluster: SE^^^ = 6E^^ + 
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X]a<b<c ^a,b,c', the noii-additivity corrections to the second-order approximation are 
ea,b,c = EQ°"{a, b, c) - [da + + Cc) - (ea,fe + ea,c + efe,c) - Eq. Continuing this way 
one builds order by order a converging series 

5Eo™- = ^e, + 5^e,,6+ 5^ ea,5,c+ ^aAc.d + • ■ ■ (3.17) 

a a<b a<b<c a<b<c<d 

which in its hmit gives the exact value of 5Eq°^^ . However, as was discussed in a num- 
ber of papers (e.g. in Ref. [Et92i], |PFS95j . |PFS96j . |KPFS97j . |DDFS97j 



and jF02j ) one obtains a very accurate estimate of 5Eq°" when the series fl3.17|l is 
truncated after the third term. Thus, in each separate correlation calculation at 
most six electrons are treated simultaneously. Moreover, as correlation effect de- 
cays rapidly with the distance between localized orbitals, the two- and three-bond 
increments and ea^^c involving distant bonds can be neglected and therefore the 
effort for evaluating the correlation correction to the ground-state energy of a big 
molecules is substantially reduced. In fact, in the case of a periodic system only a 
finite number of energy increments has to be considered. 

An extension of this approach to the correlation corrections to local SCF matrix 
elements of hole states was proposed in the paper j(TSF93j and developed further 
GSF97j . |AFS00j and |AF02j . It was shown that the corrections to the matrix 



m 



elements given by Eq. ()3.15p can be also obtained with any desired accuracy 

by means of local increments. In Eq. ()3.15p the term SH^^" accounts for correlation 
in the (TV — l)-electron system where electrons react on the presence of the hole while 
^^corr accounts for electron correlation in the neutral system. As the difference 
of these two quantities define SIP^^" they must be evaluated consistently in the 
incremental scheme. 

By analogy to Eq. ()3.17j) for the correlation correction to the ground-state energy 
the incremental approach to 51?^™"^ can be written as 



SIPT = AlPr + E ^IPr (c) + E AlPTic, d) + ... (3.18) 

cj^a,b c<d^a,b 

where in parentheses we denote those additional bonds which are "open" for corre- 
lation together with the matrix element defining bonds a and b. Below we explain 
the performance of the approach separately for diagonal and off-diagonal elements. 

The first order of approximation to the diagonal element SIP'^J^ is obtained when 
only bond a is open for correlation. This means that in the (A^ — l)-electron system 
one electron in bond a is destroyed and the second electron from this bond can 
be excited to the virtual orbitals thus forming singly excited configurations in the 
expansion (|2.76|) of the reference wavefunction: 

= «a(a)$. + J2^J,{a)^J, (3.19) 
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with energy E^°^^. In this case the effective Hamiltonian is a 1 x 1 matrix whose 
element is equal to E^"". The corresponding energy of the neutral (i.e. A^-electron) 
system is i?Q°''''(a). The first-order approximation to SIP^^^ thus is 



AlPr = SIKTO = HfA) - Haa - [Er{a) - E,] = - IPaa. (3.20) 

The second-order approximation is obtained when the contributions to (5IP^°'^'^ 
from each of the remaining (riocc ~ 1) bonds individually are summed up, where the 
corrections are defined relative to AIP™"(), 

Aipr (c) = 5ipr (c) - Aipr (3.21) 

The first term in the right-hand side of Eq. ()3.2H) is 



5IPr (c) = Hf^{c) - H^a - [Eri^. c) - E,] = IPr (c) - IPaa. (3.22) 

In Eq. (j3.22|) the effective Hamiltonian element H'^[c) stems from a 2 x 2 matrix 
defined in the configurations $a and $c- At this point we would like to notice that 
due to the decay of the correlation effect with the distance between the correlated 
electrons one finds a rapid decrease of the increments AIP^ (c) when the distance 
between the orbitals a and c increases. This ensures that one can indeed use finite 
clusters to evaluate the correlation corrections to the local matrix elements defined 
in a crystal. 

It is obvious that summing over all second-order increments does not yet give 
a correct estimate of ^IP^ as correlation between electrons in different additional 
bonds (c and d) is missing. The next order of approximation is obtained when we 
account for the corresponding non-additivity corrections 



AlPr (c, d) = 5IPr (c, d) - AlPr - [AlPr (c) + AlPr (rf)] • (3.23) 

Here the correlation correction ^IP^ (c, d) to the diagonal matrix element IPa^ is 
calculated with three bonds (a, c and d) being open: 



5IPr (c, d) = Hf^ic d) - H^, - [E™-(a, c, d) - E,] = IP-" (c, d) - IP.,. (3.24) 

The hole is delocalized within these three bonds and thus three correlated (A^ — 1)- 
electron states are obtained. This procedure can be continued further order by order. 
However, as increments are expected to decrease rapidly with increasing order (see 
Ref. |(;SF97j . |AFSnnj and |AFn2p . the series (HrTT^ can be truncated after the third 
term. 
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The correlation corrections to off-diagonal (or hopping) matrix elements IP^^" 
are evaluated the same way, however, the procedure is started with the two bond a 
and b and no ground-state energies are present: 

AlPrO = 5IPa6() = H:^0 - H^, = IFTO - IP,,. (3.25) 

Here the starting effective Hamiltonian matrix elements H^f[) already originate 
from a 2 X 2 matrix defined in configurations and Next we compute the 
correlation contribution of an individual additional orbital: 

AlPr (c) = 5lPTic) - AlPTO (3.26) 

where 

SlPTic) = H:^{c) - H,, = IPTic) - IPa,. (3.27) 

In this case the hole is delocalized over a triple of bonds (a, b and c). Summing over 
c we get the second term in the series ()3.18|) . The third-order increments are 



AlPr (c, d) = SlPTic, d) - AlPTO - [AlPr (c) + AlPTid)] (3.28) 

where correlated matrix elements 5\P'^°"{c,d) are obtained when four bonds are 
open for correlations always including the bonds a and b: 

5IPr (c, d) = H:I{c, d) - H,, = IPTic, d) - IP.,. (3.29) 

Again, the series ()3.18p can be truncated after the third term. 

In the present work this approach is used to obtain good estimates for the corre- 
lation corrections of the local matrix elements IPab- For the a bonds a CeHg molecule 
is used. Due to the truncation of the series ()3.18|) in each CI calculation we have at 
most four bonds open for correlation which keeps the calculations reasonably cheap. 
However, even in such a small molecule as CgHg the number of all third-order incre- 
ments is 78 (all pair combinations among 13 bonds) for the diagonal elements \Paa 
and 66 for the off-diagonal elements (2 out of 12 combinations). Dozens of input 
files have to be prepared for the correlation calculations and the results of these 
calculations have to be gathered and summed up properly. As discussed in the sec- 
tion 3.3 we need much larger molecules to account for the correlation corrections to 
local matrix elements \Pab with a and b being vr bonds and also for selected EA^s 
matrix elements. The number of increments increases quadratically with the size 
of the cluster (i.e. with the number of bonds). Therefore, the incremental scheme 
introduced above in terms of bond contributions becomes unfeasible because of the 
huge number of increments though each increment individually can be evaluated by 
the CI method at low cost. 
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3.4.2 Formulation in terms of bond groups 

To proceed in the case of big clusters we propose to reformulate the incremental 
scheme in terms of groups of orbitals. In this approach to the correlation corrections 
of the matrix element IPab we subdivide all bonds except of a and b in groups gi, 
g2, ... consisting of a couple of bonds each {ug. > 1, n^. being the number of bonds 
in the i-th group g^). It is natural that the bonds in a particular group gi belong to 
the same unit cell (or two adjacent unit cells) and that the groups do not overlap. 
Then, in Eq. (j3.18|) we replace the single orbitals c and d in parentheses by entire 
orbital groups gi and gj and sum over group indices i and j: 

5IP™" = AIP-"() + J2 ^^PTi9^) + E ^^^Ti9^, 9,) + ... (3.30) 

i i<j 

In the case of a diagonal element, IPaa say, the other orbital b forms its own group and 
has to be included in the summation over the groups gi in Eq. (j3.30p . Consequently, 
we replace the bonds c and d in Eqs. fl3.21|) - ()3.24|) and ()3.26|) - ()3.29|) by groups gi 
and gj. Thus, the increments in the series (l3.3Up become 

Aipr (^7.) = slKri9^) - Aipr (3.31) 

and 



AlPr (^?. 9j) = SlPTi9^, 9,) - AlPr - [AlPr (^,) + AlPr (^,)] (3.32) 

The substitution of a group gi instead of a single bond c in (j3.2ip physically means 
that we account for all the contribution to the correlation correction 5YP^°" which 
are associated with electrons from the whole group of bonds gi instead of a single 
electron from bond c: 



5\VT{9i) = H:!igi) - Ha, - S,, [E^ia, g^i) - Eo] = IPT (^?.) - IPa.- (3.33) 

To get the matrix element H^figi) of the effective Hamiltonian, {ug. + 2 — 6ab) 
correlated hole states have to be calculated explicitly. Also for diagonal elements 
the correlation correction to the ground-state energy is obtained when all bonds of 
the group gi together with the bond a are open for correlation simultaneously and the 
rest of the bonds of the cluster are frozen. Summing up the correlation contributions 
associated with the individual groups one gets a first estimate of the response of the 
cluster to the presence of a hole. Note that this first-order estimate SlP^^fj contains 
much more correlation contributions than the sum over all single bond increments 
because intragroup correlation is fully included. This result is improved when the 
non-additivity corrections AIP^^" {gi, gj) for all pairs of groups are summed up. The 
scheme can be continued by adding increments of higher orders. However, as group 
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increments also decrease rapidly with the order, the series (j3.30p can be terminated 
after the third term to get rather accurate estimate for the correlation correction to 
the LME IP ab- The incremental scheme in terms of bond groups, as it is introduced 
in this subsection, was used in our approach to evaluate the correlation correction 
to the matrix element IP±a,pi pi from Table 13.31 In that case groups of a bonds were 
composed. The corresponding data on the increments are given in chapter 4.1. 

Grouping the bonds substantially reduces the numbers of increments in ()3.30p as 
compared to (j3.18j) however to the price of more expensive CI calculations for each 
increment in (j3.3Up . In the latter case up to {ug- + Ug. + 2) bonds can be open for 
correlation simultaneously and the same number of hole states have to be calculated 
to get the matrix element H^f{gi,gj) for a particular increment AIP^^"'(5fj, (yfj). 

3.4.3 Incremental scheme for attached-electron states 

The method of local increments can be also used for the evaluation of the correlation 
corrections to the matrix elements EA^^- In this case (2nbond + 1) electrons have 
to be correlated and the additional electron is delocalized over all antibonds of the 
chosen molecule. To reduce the computational efforts for the MRCI calculations, 
one can additionally freeze the majority of antibonds when calculating a particular 
increment forbidding electrons to occupy the frozen antibonds. This way one reduces 
i) the number of calculated (iV + l)-electron states and ii) the external space. At the 
same time only electrons from the bonds belonging to one or two groups are open 
for correlation simultaneously. (Here we need to group the bonds again, because 
large clusters are required to account for correlation effects in the (A^ + l)-electron 
systems.) In the incremental scheme the extension of the external space (i.e. the 
opening of antibonds) and the opening of the bonds must be done consistently. 
Otherwise the increment series will converge poorly. Thus, we find it reasonable to 
group the virtual localized orbitals centered on one bond together with the occupied 
orbitals from this bond (or adjacent bonds). The general formula for correlation 
correction to matrix element EA^s by means of local increments is analogous to 
Eq. dSSni): 

= AEA-" + E ^^KT{9^) + E ^^^rTi9^, 9,) + ■■■ (3-34) 

where 

AEA^Ti9^) = SEA^Ti9i) - AEA™" () (3.35) 

and 

AEAr (^., 9j) = ^^AT{9^, 9j) - AEA-" () - [AEA-" (^7.) + AEA-- (^?,)] • (3.36) 
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Here indices r and s denote two groups gr and each containing the relevant 
antibond plus corresponding bond (or few adjacent bonds). Again for the diagonal 
element, EA^r- say, the summation over the groups in Eq. (|3.34|) includes the group 
gs of the other antibond. 
The first-order increment 



AEA-"() = - [H^^O - Hrs] + 5rs [Eri9r) - Eo] = EKTi) - EArs (3.37) 

contains the effective Hamiltonian defined in configurations and $s; all occupied 
orbitals from the groups g^. and gs are open for constructing the excited configura- 
tions in the MRCI ansatz for the correlated wavefunctions ()2.78p and all antibonds 
except of r and s are forbidden for the excited electrons. For the diagonal element 
(s = r) only one group is considered in Eq. ()3.37|1 . The correlated energy of the 
A^-electron state EQ°'"{gr) is calculated with only bonds of the group gr being open. 

When we add individual groups for the second-order increments and pairs of 
groups for the third-order increment we simultaneously open antibonds and the 
corresponding bonds for correlation. This way we approach in a systematic way the 
value of (5EA™'^'^ for all bonds and all antibonds being open simultaneously, however, 
in each separate MRCI calculation the number of configurations (and the number 
of coefficients a in Eq. (|2.68j) and ()2.78|) to be optimized) is considerably restricted. 



3.4.4 Additional features of the incremental scheme 

The incremental scheme described above has five important additional features. The 
calculation of individual second-order increments A\P^°j^^{gi) and ^EA^°^^{gi) gives 
information on the decay of the correlation effect with the distance between the 
bonds of a group gi and the localized charge in bonds a and h (or antibonds r and 
s). This information allows one to estimate whether a chosen cluster is big enough 
to account for the dominant part of the correlation corrections. 

The second feature is that one uses the symmetry of the crystal to reduce the 
number of increments which have to be calculated: the contribution of symmetry- 
equivalent increments to the correlation correction is the value of one of them times 
the number of such equivalent increments. This reduces the computational efforts 
and works especially well for highly symmetric crystals with simple unit cells like 
diamond, silicon and germanium (see Ref. |(;SF93j . |(;SF97j . |AFSnnj and |AFn2j 1. 



The third feature: the main computational effort is directed to get the effective 
Hamiltonian from which we only need one of its matrix element for a particular 
increment. However, it is easy to see that some other matrix elements of a given 
effective Hamiltonian can be used for other increments without any extra cost. For 
instance, the diagonal element if°f(fe, c) and the off-diagonal element H^f{c) are 
obtained from one and the same MRCI calculation with $„, $6 and $c being the 
model configurations. 
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The fourth feature is the possibihty to use the translational symmetry of the 
crystal and reduce the so-called cluster- edge error by this. Such errors arise when 
one of the bonds open for correlation is close to the edge of a cluster and this localized 
orbital is affected by the lack of supporting atoms (i.e. the energy of the electron 
in this orbital differs noticeably from the energy of the electron in a translation- 
equivalent orbital in the central part of the cluster). Then, for the evaluation of this 
particular increment one can open a translational-equivalent set of orbital groups 
for correlation which does not contain orbitals at the edge of the cluster if the size 
of the cluster allows one to do so. Let's say that we want to get the contribution to 
the correlation correction to the diagonal matrix element IPpi,pi associated with the 
additional bond "pi,-2" in the molecule C10H12 which is shown in Fig. 13.71 We know 
that orbital "pi,-2" is not equivalent to orbital "pi,-!" as one more unit cell from the 
left side is needed to reproduce properly the occupied Wannier vr orbital centered on 
the bond "pi,-2" (see Fig. 13. 6j) . Also we expect that the increment AIPp°pj(pi, —2) 
is still relatively large as orbitals "pi" and "pi,-2" overlap as can be understood from 
Fig. 13.61 Therefore, this increment has to be evaluated when two equivalent bonds 
with proper crystalline surrounding are open for correlations and in the series ()3.18|) 
we substitute it by the translation-equivalent increment AIPp°[''pj]^(pi, —1). This way 
we reduce the influence of the edge of the cluster on the correlation corrections to 
the local matrix elements defined for the periodic infinite systems. 

The fifth feature becomes important when one uses a correlation method which 
suffers from size-consistency problems. This is particularly the case when truncated 
CI methods are used. When the second- and third-order increments are calcu- 
lated, correlation corrections obtained with different numbers of open bonds are 
subtracted from each other. The difference itself may easily be of the order of the 
size-consistency error (better to say size-extensivity error) for each term entering 
this difference. Then, by the incremental scheme, we would sum up these errors 
which would accumulate and finally lead to wrong results. Therefore, one needs 
an accurate correction to minimize the size-extensivity error when the MRCI(SD) 
method is used for the evaluation of the increments. Such a correction is the topic 
of the next section. 

3.5 Size-extensivity correction 

The occurrence of a size-extensivity error in truncated CI calculations can easily 
be understood from the example of two H2 molecules at large distance from each 
other. Each molecule has one bond (an orbital obtained in the Hartree-Fock ap- 
proximation) occupied by two electrons. Let us evaluate the correlation correction 
SEq°"{2H2) to the HF ground-state energy Eq{2H2) of this system by the single 
and double configuration interaction method (for shortness we will call here this 
correction as correlation energy). We can do this in two ways. Firstly, we can 
calculate the correlation energies for each H2 molecule separately, 6Eq°"{H2), and 
add them. This approximation gives exact result since on big distances between 
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the H2 molecules excited configurations with electrons being destroyed in the bond 
of one molecule and created in some virtual orbital of the second H2 molecule give 
no contribution to the correlation energy (dissociation limit). Alternatively, we can 
directly calculate the correlation energy of the whole system of two H2 molecules. 
Due to the truncation of the CI expansion ()2.52|) we have more configurations in the 
first case than in the second one: for instance, double excitations on both molecules 
simultaneously are neglected in the second case while they are taken into account 
in the first case. Therefore, the correlation energy obtained in the whole system 
5Eq°^^{2H2) is higher than the sum of the two correlation energies of the separate 
H2 molecule 25Eq°^^{H2). The error in the CI(SD) correlation energy of the system 
consisting of n H2 molecules increases with increasing n (see the third column in 
Table ESI). Thus, when one investigates correlation energy of systems with different 
number of bonds open for the correlation one has to correct for the size-extensivity 
error. 

Explicitly this situation arises when one uses the incremental scheme to evaluate 
correlation energies. There, correlation energies of the system with different numbers 
of correlated bonds are subtracted from each other. Therefore, all the energy values 
defining a particular increment must be consistent. 

The general strategy to obtain size-extensive correlation energy is to estimate 
the value n x 5Eq°^^{H2) by means of information on the wavefunction and the 
correlation energy 6EQ°"{nH2) of the whole system. There exist two well-known 
size-extensivity corrections to the correlation energy of the ground state calculated 
by the CI method which are widely used in quantum chemistry. 

The first is the Davidson correction |LD74j which is the simplest estimate of the 
size-extensivity error in CI(SD) calculations: 



Here ao is the coefficient of the SCF ground-state configuration $ in ()2.68|) when 
the correlated wavefunction is normalized. The Davidson correction does not 
depend explicitly on the number of correlated bonds. It gives a rather accurate 
estimate of correlation energy for n ranging from 4 to 10 (the residual error is 
within 1% as pointed out in |S096j . p. 268) however the error is larger for n = 2 
or 3 (see Table ITBj) and it diverges for large numbers of open bonds. Since in our 
incremental calculations we have dozens of increments the error may accumulate up 
to 100% of the correlation energy. Thus, for our purposes we need a more accurate 
size-extensivity correction. 

The Pople correction |PSK77j estimates the size-extensivity error much better 
than the Davidson correction (see Table giving perfect size-extensive correlation 
energy in the case of any number of equivalent non-interacting two-level systems: 



AEf^ 



D — 



5Er{l-al). 



(3.38) 




(3.39) 
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Figure 3.14: Four H2 molecules at large distances. The contour plot of the bonding 
HF orbital of one H2 molecule is shown. 

Table 3.6: Comparison of the size-extensivity errors of the CI(SD) method with 
and without error corrections. The correct size-extensive correlation energy for 
different numbers n of correlated H2 molecules is presented in the second column 
(in eV). The correlation energy of three different methods (pure Cl(SD), CI(SD) 
with Davidson correction and CI(SD) with Pople correction) are presented in the 
three next columns as deviation from the correct correlation energy. 



n 


n 5^0"" (H2) 


5Eo™"(nH2) 


6E^°''{nR2) + A^SE 


6E^°''{nR2) + ^Ef^ 


2 


-2.063 


1.721% 


-1.706% 


0.008% 


3 


-3.094 


3.318% 


-1.533% 


0.031% 


4 


-4.126 


4.806% 


-1.318% 


0.064% 



where 6 = arccosao. 

We have investigated the reduction of the size-extensivity error of the CI(SD) 
method in evaluating the correlation correction to the ground-state energy when the 
two error corrections are used. As a test system we have chosen four H2 molecules in 
the arrangement shown in Fig. 13.141 In the dissociation limit the distance between 
nearest-neighbor molecules (20 A) is much larger than the distance between the two 
H atoms in a molecule (the characteristic bond length in H2 molecule is 0.7414 A 
which is taken from Ref. |HW76j . p. 52). The basis set for hydrogen used here is 
the same as in Table ITTl The correlation energy of the system with n correlated H2 
molecules (n = 2, 3 and 4) obtained by the CI(SD) method is compared in Table ITHl 
to the size-extensive one which is equal to n times the correlation energy of a single 
H2 molecule. 

The Pople correction reduces the size-extensivity error in average by two order 
of magnitudes as compared to pure CI(SD) or Davidson corrected CI(SD) which 
both yield error of more than 1% (see Table ESI). Thus, we use the Pople correction 
to Cl(SD) correlation energies to get size-extensive results for the neutral closed- 
shell systems. (More detailed study of the performance of these two corrections for 
different systems can be found in e.g. Ref. |PSK77j . |LBT87i . |C;A88j and IDDMj ). 
Unfortunately, being applied to the correlation energies of open-shell systems (as the 
{N — 1)- and (A^ -|- l)-particle systems we need to consider) the two size-extensivity 
corrections do not yield satisfactorily accurate results (see Table 13. 7j) to be used in 
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the method of local increments. Therefore, we need a size-extensivity error correc- 
tion which gives accurate estimates of the correlation energy in the case of open-shell 
systems (systems with a hole or with captured electron) like the Pople correction 
does for closed-shell systems. 

Let us consider again the system consisting of n equivalent H2 molecules at large 
distances and let us consider a {2n — l)-electron state with one electron annihilated 
at some molecule a. We can get the correlation energy of this state on the MRCI(SD) 
level but we would like to approach the correct size-extensive correlation energy 

SE^ = (n - 1) SErm + SEnnt) (3-40) 

having information on the correlation energy and the wavefunction from the MRCI(SD) 
calculation of the whole system. In Eq. ()3.40j) 6Eq°'^'^{B.2) denote the correlation en- 
ergy of a neutral H2 molecule and 5£'™''''(H^) the one of the molecule carrying the 
hole. The latter is in fact just the orbital relaxation energy. 

Like in the derivation of the Pople correction |PSK77j we switch to a minimal 



basis description such that the H2 molecule becomes a two-level system. The spin- 
adapted reference configurations for the MRCI ansatz for our system consists of 
all hole configurations {^a}a=i,...,n with a = 1, . . . ,n being the occupied spin-down 
orbitals on the different H2 molecules. In the dissociation limit the only 2-hole 
1-particle configurations which couple to a given model configuration <I>a, i.e. 
for which the CI matrix element is non-zero, are the "vertical" single 

excitations $^^* where b* labels the unoccupied counter part of the spin orbital ipb. 
Because of Brillouin's theorem ) vanishes if the spin orbital ipb is located 

on a different molecule than a. Hence, only the configuration $^"* couples to $a. 

Similarly, the only 3-hole 2-particle configurations ^^cd which couple to a given 
model configuration are the "vertical" double excitations ^J'^i* on sites dif- 
ferent than the one of the spin orbital (fa- Furthermore, the configuration space 
{'^a}a=l,...,n® {^aa }a=l,...,n © {^J'bf}<^^f>=i,...,n is closed in the sense that no other 2- 
hole 1-particle or 3-hole 2-particle configuration couple to any of the configurations 
in that space. Hence, (in the dissociation limit) we can write the MRCI wavefunction 
of our system as 



a=l a=l a=^b=l 

For the sake of simplicity let us denote $„ by <I>°, $J* by ^l, ^J"^^* by . . .. 
Then, the CI matrix elements for the MRCI ansatz ()3.41|) can be written as 

c:i: = mH\<^i). (3.42) 

In the dissociation limit the CI matrix becomes block-diagonal, i.e. 
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with all subblocks {C^a}x,y=o,...,n being of the form 
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(3.43) 



The entries Eih, E2h, E^h, P and 7 do not depend on the position of the hole in 
the H2 chain which implies that our {2n — l)-electron system of nB.2 molecule has 
n degenerate lowest hole states each of them residing on one of the H2 molecules. 
Delocalization of the hole over the chain does not lead to a gain in energy. 

The matrix elements 7 = {^alHl^J'*^ ) = (ipb'^blvul'^l'^l) are identical for all 
b ^ a since all the neutral H2 molecules are equivalent (in the dissociation limit). 



For that reason the CI coefficients a 



b*b* , 
a bh 



must all be equal, and the wavefunction 



for a hole state on the a-th H2 molecule (with energy E'^^^) can be written as 



^corr 
a 



(3.44) 



x=2 



As a consequence the secular equation to determine the coefficients Ci and C2 reduces 
to 
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(3.45) 



where SE''^" 



_ ^corr _ ^YiQ correlation energy, and 2Ai = E2h — Ea and 

2A2 = E'ih — Ea are the energies of the 2-hole 1-particle and the 3-hole 2-particle 
configurations (^^n ^"^^ ^all ' respectively, relative to the Hartree-Fock energy 
Eih = Ea of the model configuration $a. From the first row one concludes that the 
correlation energy consists of two parts: SE^"^^ = ci(3 + (n — l)c27 = Es + Ed where 
the part Es (Ed) is associated with singly (doubly) excited configurations only. The 
quantities Eg and Ed (and also = Ci and q = y/n — 1 C2) can be obtained from the 
output of an MRCI(SD) calculation with the MOLPRO program packadge |M2000j . 
|WK88j . |KW88j . |KW92j and our aim is to express the size-extensive correlation 
energy in terms of these parameters. 

As can be seen from Eq. (I3.40|) the size-extensive correlation energy of the con- 
sidered system is equal to the correlation energy the of charged molecule plus 
(n — 1) times the correlation energy of a separated neutral H2 molecule. To find the 
required energies 5i?™''''(H^) and 6Eq°"(B.2) of the two isolated molecules (charged 
and neutral) we note, that one gets the following CI matrices when the isolated 
molecules are treated analogously to the whole system: 



(3 
P 2Ai 



1 



sEn^t 



(3.46) 
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and 



7 
7 2A2 



5E',"''{R2) 



(3.47) 



with the entries /3, Ai, 7 and A2 being the same as in Eq. ()3.45j) . From these two 
equations the correlation energies of the charged and neutral H2 molecules can easily 
be derived yielding: 



and 



Ai 



A, 



A^ + y. 



(3.48) 



(3.49) 



Substitution of ()3.48j) and ()3.49j) into Eq. ()3.4()|1 gives the size-extensive correlation 
energy of the lowest n hole states in our system of n H2 molecules in terms of the 
parameters entering the secular equation ()3.45|) : 



<5i?f = Ai - y/A? + /52 + (n-l)(^A2 - ^ + 7^^ . 



(3.50) 



With the help of Eq. ()3.45j) the parameters P, 7, Ai and A2 can be linked to the 
output quantities SE^°", Es, Ed, Cs and q which are available from the MRCI(SD) 
calculations for the n x H2 molecules whole: 



2Ai 
2A2 



6Er- 



Es/cl 



7 = Ed/V^^ Cs , 2A2 = 5^^=^ - Ed/ci . (3.51) 
Then, the final formula for the size-extensive energy of the hole state reads 



SE': 



■BE 



1 

2 

(n 



n5Er-^-{n-l)^ 
ci ci 





Ed 
^2 



Ej 



{n - l)cl 



(3.52) 



This correlation energy effectively takes into account the missing configurations in 
the truncated CI expansion of the wavefunction of an open-shell system as a function 
of the number of bonds open for correlation. For a system more complicated than a 
chain of isolated two-level systems the quantities Eg (Ed) and (q) are associated 
with the contribution from all single (double) excitations to the correlation energy 
E^"'^'^ and the squared norm of the (normalized) wavefunction respectively, as 
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Table 3.7: Comparison of the size-extensivity errors of the CI(SD) method with and 
without error correction (three different methods) for a hole states. The correct 
size-extensive correlation energy for different numbers n of correlated H2 molecules 
is presented in the second column (in eV). The relative deviation of the correlation 
energy of four different methods (the same as in Table 13.61 plus our own correction 
SE^^) are presented in the four next columns. 



n 




5E™-(nH2) 


6E^^°''{nE2) + AEl^ 


5E^°"(nH2) + AE|^ 


6E!^ 


2 


-1.857 


2.342% 


-2.176% 


2.342% 


0.002% 


3 


-2.888 


4.050% 


-1.749% 


1.153% 


0.017% 


4 


-3.919 


5.553% 


-1.400% 


0.807% 


0.045% 



is also done usually for the Davidson and Pople corrections to the MRCI correlation 
energies in quantum chemistry. 

To test this new size-extensivity correction for the MRCI(SD) method we per- 
form correlation calculations for the hole states in our 4 x H2 system. The correlation 
energies of this system size-corrected by different methods (Davidson and Pople cor- 
rection as realized in MOLPRO jM2000j . jWK88j . |KW88j . jKW92j and 6E^^ from 
Eq. ()3.52|l ) are compared to the correct size-extensive correlation energy computed 
for a separate charged molecule and (n — 1) neutral H2 molecules (denoted as 
^^exact shortness). They are presented in Table ITTI for different number of open 
bonds (n = 2, 3 and 4). The distance between nearest-neighbor H2 molecules in 
these calculations is set to 50 A to reach the dissociation limit in our charged system. 

From Table ITTI one concludes that the new size-extensivity correction for open- 
shell systems proposed in this section gives results which are by several orders of mag- 
nitude more accurate than the convensional corrections to the MRCI(SD) method. 
The accuracy obtained with our size-extensivity correction is the same as that of the 
Pople correction to correlation energy of a neutral closed-shell system. The use of 
these two corrections in all our calculations for individual increments (the first one 
for open-shell and the second for closed-shell calculations) finally allowed us to re- 
duce the estimated error for the correlation corrections to the local matrix elements 
(the target quantities of our approach) down to a few percents of their values. 

3.6 Satellite states problem 

In the last section of this chapter we would like to draw our attantion to the problem 
of satellite states which may ruin the convergence of the MRCI calculations. This 
problem appears when one (or more) of model configurations $j which define the 
effective Hamiltonian H^^ for some particular increment has an energy higher than 
the energy Ejjl of a single excitation of another model configuration $j. This 
situation may happen in systems with relatively broad bands and a narrow band 
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LUMO 7t 

71 hole a hole 

Figure 3.15: Three possible configurations for an (A^ — l)-electron states of tPA: 
single hole in iether a vr or a a bonds and a satellite state of vr symmetry. 

gap. For instance, in the case of trans-polyacetylene the HF energies of hole states 
range from approximately -31 to -6 eV and the band gap is 6.42 eV. The position 
and dispersion of the lowest valence band of tPA (see Fig. 13. 4|) is mainly defined by 
states with a hole in the carbon sceleton. In a CgHg molecule the lowest localized 
valence orbital is a short C-C a bond and the energy of taking out an electron from 
this bond is 22.3 eV. The position and the dispersion of the highest valence band 
(a vr band) in turn is determined by states with a hole in the vr system of tPA. In 
the CgHg molecule the highest localized valence orbital is a vr bond with an orbital 
energy of 10.4 eV. The energy of a state with two electrons destroyed in that vr 
bond and one electron created in the lowest unoccupied molecular orbitals (LUMO) 
is found to be 19.8 eV in a CeHg molecule, and is lower than the energy of the a 
hole. A schematic sketch of this three configurations is shown in Fig. 13.151 

If in such a case we would perform MRCI(SD) calculations to determine the 
correlation contributions to a matrix element IP„„{7t), say, we would construct the 
model space of two one-particle configurations and $o- and calculate the two 
lowest correlated (A^ — l)-electron states with wavefunctions defined by Eq. ()2.76|1 . 
During the CI iteration a decrease of the CI coefficient acr{2) for the model configu- 
ration $o- and an increase of the CI coefficient a^'^* (2) for the satellite configuration 
(which is present in the CI expansion ()2.76|) ) would lead to a gain in energy and 
after few iterations we would arrive at a state which is dominated by the satellite 
configuration $^^* instead of the model configuration The obtained correlated 
wavefunction does not obey the requirement of the MRCI approximation that ref- 
erence configurations (the model configurations here) must be dominant in the CI 
expansion of the wavefunctions. At this point the program stops the iteration with 
an error message and no data for the two correlated states are obtained. To pro- 
ceed in this situation we include the satellite configuration in the reference space 
and calculate the first three correlated states. Two hole and one satellite state will 
emerge from this calculation. 

The hole states can easily be separated from the satellite state by checking the 
square of the norms {aa{i)Y of their projection onto the subspace spanned by the 
one-particle model configurations. In all our calculations where satellite states show 



H 

% satellite 



E'^'< E 



74 



CHAPTER 3. REALIZATION OF THE METHOD 



Table 3.8: Correlation data for two hole and one satellite state in a CqHs niolecule. 
Correlation energies are given with respect to the HF ground-state energy of the 
molecule. 



state 


projected 


E^"" (eV) 


71 hole 


0.96 


8.63 


TT satellite 


0.00 


16.15 


a hole 


0.86 


18.54 



up the squares of the norms of the projected a hole states are higher than 0.70 such 
that the distinction between satellite states and single hole states is clear. Thus, the 
requirement that we have to find those wavefunctions ^1'^°'''' which are dominated by 
model configurations is fulfilled. 

In Table EIHl we present the correlation data corresponding to the case described 
above. In the second column the squares of the norm of the correlated wavefunctions 
projected onto the subspace spanned by model configurations <I>^ and are shown. 
The wavefunction of the satellite state has essentially zero projection and is clearly 
more stable than the deep-lying a hole state. 

In the present case, the inclusion of singly excited configurations into the ref- 
erence space helps to solve the problem of competing states in MRCI calculations 
and allows us to get the desired target matrix elements of the effective Hamiltonian. 
However, this approach has limitations. When many one-particle configurations 
with relatively large energy dispersion are included in the model space or bigger 
clusters are considered (which implies more low-lying unoccupied orbitals) there 
may be a lot of satellite states competing energetically with the one-particle states. 
Therefore, many states have to be optimized simultaneously in the MRCI iterations 
and such calculations may become quite expensive (or even unfeasible). To pro- 
ceed in this case one needs an improved version of the MRCI program which could 
sort out the relevant states at each step of the iteration (again by the norm of the 
projection) and converge only states with a dominant contribution of one-particle 
configurations in the CI wavefunctions. Such a treatment would completely solve 
the problem of satellite states in the MRCI method as long as the satellites do not 
strongly mix with the single hole states. 



Chapter 4 

Results for ^rans— polyacetylene 



In this chapter we present and discuss the numerical results for trans-polyacetylene. 
In section 4.1 the individual increments obtained in clusters are given and the con- 
vergence and accuracy of the incremental scheme is explored. The correlation cor- 
rections to the local matrix elements are obtained and the change of the LMEs 
due to electron correlation is discussed. In the next section we present the corre- 
lated band structure of trans-polyacetylene and compare it to that on the HF level. 
In section 4.3 we estimate the error for the individual local matrix elements, the 
band gap and the band widths. In section 4.4 we discuss how results obtained for 
the one-dimensional system can be compared to experimental data on bulk trans- 
polyacetylene. 



4.1 From single increments to final local matrix 
elements 

We start the discussion of the numerical results with the correlation corrections 
to the local matrix elements associated with a hole in a bonds. In this case, the 
convergence of the correlation effect with the distance from the localized hole is 
most rapid as a bonds are the most compact localized orbitals. Therefore, we only 
need a relatively small CeHg cluster to account for the correlation effects. 

First, we would like to discuss the effect of truncating the incremental scheme 




Figure 4.1: CeHg cluster of tPA with bond labels. 
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after the third order. For this purpose we perform a correlation calculation for 
nine hole states when the most inner nine bonds of CgHg are open simultaneously 



("CH1,-1", "CC", "CH2", "sig", "CHI", "CC,1", "CH2,1", "sig,l" and "CHI,!"). 



The names of the bonds are indicated in Fig. 14.11 We do not include tt bonds here to 
avoid the appearance of satellite states. From this calculation we get the correlation 
corrections to the local matrix elements which are associated with these nine bonds 



We can get the same corrections by the method of local increments. The cor- 
responding data on the individual increments are presented in Tables 14.11 and 14.21 
where increments of different orders are placed in different columns. Numbers at the 
bottom of these tables are the sums over all increments AIP^'^'' of the given order 
n. These numbers demonstrate the convergence of the incremental scheme with the 
order of increments: the third term in Eq. ()3.18p drops with respect to the second 
one by more than an order of magnitude that insures the possibility to truncate 
the series (I3.18|) after the third term. Actually, non of the third-order increments 
of the off- diagonal metrix elements is larger than 5 meV, so that these contribu- 
tions might be neglected as well. Summing up the increments of 1st, 2nd and 
3rd order one gets the correlation correction estimated by the incremental scheme: 
'^"IPcHi.CHi = -1-149 eV and 5"IPchi,cc,i = -0-389 eV. The differences between 
the exact and the estimated values of correlation corrections are —72 meV (or 6.3%) 
and —18 meV (or 4.6%) for the diagonal and off-diagonal matrix elements, respec- 
tively. Thus, we can argue that the incremental scheme estimates the correlation 
corrections with an error of a few percent. 

From Tables 14.11 and 14.21 one also nicely sees the convergence of the second- 
order increments with the distance to the added a bond c. To study this decay of 
correlation effect with the distance from the localized extra charge more extensively 
we provide the second-order increments of the matrix elements IPcc.cc (diagonal) 
and IPccsig (off-diagonal) in Table 14.31 ordered by the distance of the added bond 
from the reference bonds. The largest value occur for the added bond adjacent 
to the one of reference bonds (or even to both in the case of off-diagonal matrix 
elements). Values of increments drop by more than an order of magnitude when the 
added bond is moved outwards by one unit cell. The contribution to the correlation 
correction coming from farther bonds is negligibly small and will not be treated 
explicitly by our correlation calculations. Rather the infinite-size extrapolation of 
the correlation correction obtained in a cluster can be estimated in a continuum 
approximation (see section 4.3). The data presented in Table justifies the use 
of the CgHg molecule as the smallest cluster of irans-polyacetylene to calculate the 
correlation corrections to all local matrix elements corresponding to cr hole states. 

Of course, the vr orbitals also contribute to the correlation corrections of the local 
matrix elements of cr-type. Together with the pure a- and mixed Tr-cr-type contri- 
butions they sum up to the total correlation corrections summarized in Table 14.41 
where the results for all diagonal and first-nearest-neighbor off-diagonal matrix ele- 




corr 



CHl.CHl ~ 



1.077 eV and one off-diagonal element 
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Table 4.1: Increments for the correlation correction to the diagonal matrix element 
IPcHi,CHi of a CgHg molecule (in eV) when only the a bonds are open for correla- 
tion. The reference value S'IV^q^^ where the electrons from all nine a bonds are 
correlated is given as well. 



^IPcHl.CHlO ^IPcHl,CHl('^) ^IPcHl,CHl(c7 C^) 



c 


AlP 


c 


d 


AlP 


0.188 sig 


-0.336 


sig 


CC,1 


0.013 


CC,1 


-0.393 


CH2 


sig 


-0.064 


CH1,-1 


-0.069 


CC 


sig 


-0.012 


CH1,1 


-0.073 


CH1,-1 


sig 


-0.000 


CH2 


-0.123 


CC,1 


CH2,1 


-0.070 


CH2,1 


-0.114 


CC,1 


sig,l 


-0.015 


CC 


-0.092 


CC,1 


CH1,1 


0.003 


sig,l 


-0.080 


CH2 


CC,1 


0.013 






CC 


CC,1 


0.008 






CH1,-1 


CC,1 


0.009 






CH1,-1 


CC 


0.003 






CH1,-1 


CH2 


0.003 






CH1,-1 


CH2,1 


0.006 






CH1,-1 


sig,l 


0.003 






CH1,-1 


CH1,1 


0.005 






CC 


CH2 


-0.017 






CC 


CH2,1 


0.007 






CC 


sig,l 


0.003 






CC 


CH1,1 


0.005 






CH2 


CH2,1 


0.011 






CH2 


"1 

Slg,l 


0.006 






CH2 


CH1,1 


0.006 






sig 


sig,l 


0.005 






sig 


CH2,1 


0.008 






sig 


CH1,1 


0.006 






CH2,1 


sig,l 


-0.010 






CH2,1 


CH1,1 


0.002 






sig,l 


CH1,1 


0.007 


0.188 


-1.280 






-0.057 


y^3 ATp(«) 
Z^i=l ^■'■■'^ CHI, CHI 








-1.149 


■'■-'^ CHI, CHI 








-1.077 
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Table 4.2: Increments for the correlation correction to the off-diagonal matrix el- 
ement IPcHi,cc,i of a CeHg molecule (in eV) when only the a bonds are open for 
correlation. The reference value S'IPq^i qq i where the electrons from all nine a 
bonds are correlated is given as well. 



^■'■PcHl.CC.lO ^-'^PcHl,CC,l('^) ^-'^PcHl,CC,l('^' 



c 


AIP 


c 


d 


AIP 


-0.213 sig 


-0.118 


sig 


CH1,-1 


0.001 


CH2,1 


-0.029 


sig 


CC 


0.004 


sig,l 


-0.024 


sig 


CH2 


0.004 


CH1,1 


-0.014 


sig 


CH2,1 


-0.001 


CH2 


-0.002 


sig 


sig,l 


0.001 


CH1,-1 


0.001 


sig 


CH1,1 


0.004 


CC 


0.000 


CH1,-1 


CC 


-0.002 






CH1,-1 


CH2 


0.000 






CH1,-1 


CH2,1 


0.000 






CH1,-1 


sig,l 


0.003 






CH1,-1 


CH1,1 


0.001 






CC 


CH2 


0.000 






cc 


CH2,1 


0.003 






CC 


sig,l 


0.001 






CC 


CH1,1 


0.001 






CH2 


CH2,1 


0.002 






CH2 


sig,l 


0.002 






CH2 


CH1,1 


0.000 






CH2,1 


sig,l 


-0.004 






CH2,1 


CH1,1 


-0.003 






sig,l 


CH1,1 


-0.005 


-0.213 


-0.186 






0.010 



-0.389 
-0.407 
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Table 4.3: Convergence of the second-order increments of the matrix elements 
IPcc,cc and IPccsig with the distance to the added a bond (in eV). 



c 


AIPcc';cc(c) 


c 


Aips'(5:.ig(c) 


CH2 


-0.576 


CH2 


-0.129 


sig 


-0.380 


CHI 


-0.064 


CHI 


-0.132 


CC,1 


-0.031 


CH2,1 


-0.025 


CH2,1 


-0.007 


sig,l 


-0.022 


sig,l 


-0.003 



Table 4.4: Correlation corrections to local matrix elements corresponding to a bonds 
(in eV). The last column gives the total sum of the incremental contributions. 



matrix element 


AIP^O 


E,AIP^-(c) 


Ec<.AiP:;^^(c,rf) 


5IP^-^ 


IPcHl,CHl 


0.188 


-3.274 


0.468 


-2.804 


IPcc,cc 


0.166 


-5.422 


1.554 


-3.702 


IP ■ ■ 


0.064 


-3.834 


-0.098 


-3.867 


IPcHl,CC,l 


-0.213 


-0.514 


0.011 


-0.718 


IPcHl,sig 


-0.210 


-0.816 


-0.014 


-1.040 


IPcCsig 


-0.165 


-0.756 


-0.014 


-0.936 



ments are given. They are obtained by the incremental scheme when electrons from 
all bonds of the CgHg molecule are correlated: the a bonds discussed above but 
also the vr bonds. Only the bonds with saturating hydrogen atoms are excluded. 
To show the convergence of the incremental scheme we provide the data for each 
increment order of the series ()3.18|) in different columns of this table. As a rule, the 
major contributions are coming from the second-order increments AIP^^''''(c). The 
third-order increments AIP'^"{c,d) are usually more than an order of magnitude 
smaller than the second-order ones, and even in the worst case (IPcc,cc) they are 
by a factor of 4 smaller. 

The same cluster is used to estimate the contribution from the a bonds to the 
correlation correction of the diagonal matrix element IPp"";. However, one needs 
CsHio molecule to get useful values for the correlation correction to the off-diagonal 
matrix element IPpi'jJi i as the two tt bonds of trans-polyacetylene are spread over 
four unit cells. In this case the adjacent cr bonds were grouped in pairs as indicated 
in Fig. O 

Nevertheless even this cluster turned out to be not big enough for the evaluation 
of the contributions from the vr bonds of tPA to ^IPpj^jJ; and (JIPpj^'p; ^. The localized vr 
orbitals are more diffuse than the a ones (see Fig. 13. 6|) and therefore the correlation 
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Figure 4.2: Grouping of a bonds in the CgHio cluster of tPA for the calculation of 
the correlation correction of IPpi,pi,i matrix element. 

effects decay noticeably slowler with the distance from the it hole. Hence, the 
calculation of the increments with added vr bonds is performed in a C10H12 cluster 
which contains five bonds of tt symmetry (see Fig l3.7p . 

We are going to show again the accuracy of the incremental scheme in the case 
of vr bonds by first calculating the correlation corrections to the matrix elements 
IPpj^jJi and IPpi^^i 1 when all five vr bonds of C10H12 are correlated simultaneously 
^ _i;i45 eV and 6'lPl°^^-^^^ = -0.573 eV ) and comparing them to the 
values obtained by the incremental scheme. The data on the increments are pre- 
sented in Tables 14.51 and 14.61 The convergence of the incremental scheme with 
the order is clearly discernable. Yet, there are also some third-order contributions 
which are still in the order of the second-order terms (~0.3 eV). Those are not 
increments AIPpi^pi(c, c?) with the orbitals c and d being arranged symmetrically 
around the reference vr bond (as one would intuitively expect) but a linear config- 
uration of the three bonds with the two added bonds c and d on the same side, 
i.e. AlPpi pi(pi, 1; pi, 2). Summing up the individual contributions in each table one 
gets the total correlation corrections for the five 7r-bonds as estimated by the incre- 
mental scheme: S"!?""";^-^ = -1.358 eV and ^'TP^-p^j ^ = -0.615 eV. The differences 
between the exact and the estimated values of correlation corrections are 213 meV 
(or 15.7%) and 42 meV (or 6.8%) for the diagonal and off-diagonal matrix elements, 
respectively. We see that the performance of the incremental scheme is still good 
even in the case of rather diffuse orbitals and in the worst case it produces an error 
of approximately 0.2 eV. 

To illustrate the decay of the correlation effects with the distance from the local- 
ized vr hole a representavive selection of the second-order increments obtained in the 
C10H12 molecule when only vr bonds are added are given in Table 14.71 The values 
of the increments of the diagonal matrix element drop by an order of magnitude 
when distance of the added bond increases by two lattice constants. The decay of 
the increments for the off-diagonal matrix element is even more rapid. The dash in 
Table l^n f and also in Table means that the absolute value of the corresponding 
increment is below the estimated size-extensivity error for this increment. 

The increments for the correlation corrections to the matrix elements IPp""; 
and IPpi pj ]^ obtained on the CeHg and C10H12 clusters are summed up for both. 
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Table 4.5: Increments for the correlation correction to the diagonal matrix element 
IPpi.pi of a C10H12 molecule (in eV) when only the vr bonds are open for correlation. 
The reference value where the electrons from all five vr bonds are correlated 

is given as well. 



^jpcorr Q 
pi,piV 



^jpcorr 
pijPi 



AlP 



AlP--(c,d) 



d 



AIP 



0.445 



pi,l 
pi,2 
pi,-l 
pi,-2 



-0.604 
-0.162 
-0.604 
-0.162 



pi,2 

pi,-2 

pi,-2 

pi,2 

pi,2 

pi,l 



pi,l 

pi,-l 

pi,l 

pi,-l 

pi,-2 

pi,-l 



-0.291 
-0.291 
0.061 
0.061 
0.036 
0.154 



0.445 



-1.532 



-0.271 



=1 '—^^^ pi, pi 



r'jpcorr 

pi, pi 



-1.358 
-1.145 



Table 4.6: Individual increments for the correlation correction to the off-diagonal 
matrix element IPpi,pi,i of a C10H12 molecule (in eV) obtained when only the vr bonds 
are open for correlation. 



aip™:p\i() 


AIP 


corr ( „\ 
pi,pi,lvW 


AlPplpi.i 


[c, d) 




c 


AIP 


c 


d 


AIP 


-0.460 


pi,2 


-0.082 


pi,-l 


pi,-2 


-0.043 




pi,-l 


-0.079 


pi,-l 


pi,2 


0.031 




pi,-2 


-0.005 


pi,-2 


pi,2 


0.015 


-0.460 




-0.157 






0.003 


l^i=l ^^^pi,pi,l 








-0.614 


r/Tpcorr 
" ■'-'^pi,pi,l 










-0.573 
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Table 4.7: Convergence of the diagonal (left) and first nearest-neighbor off-diagonal 
increments AIPl°^^{c) with the distance between the additionally correlated bond 
(in eV). The dash means that the increment is too small to be reliable due to size- 
extensivity problem of MRCI. 



a, a 


c 


AIP 


a, b 


c 


AIP 


pi,pi 


pi,-l 


-0.604 


pi,pi,l 


pi,-l 


-0.079 


pi,pi 


pi,-2 


-0.162 


pi,pi,l 


pi,-2 


-0.005 


pi,l,pi,l 


pi,-2 


-0.053 


pi,l,pi,2 


pi,-2 




pi,2,pi,2 


pi,-2 


-0.015 









Table 4.8: Correlation corrections to the local matrix elements IPp°pi and IP' 
(in eV). The total values are presented in the last column. 



matrix element 










IP ■ 

pi, pi 


0.445 


-2.006 


-0.446 


-2.007 


TP 

-LJ- pi,pi,l 


-0.460 


0.087 


0.057 


-0.316 



cr and vr orbitals, and presented in Table 14.81 The contribution of the part of the 
crystal outside the C10H12 cluster will be estimated in section 4.3 using the data 
from Table 14.71 There we show that even in the case of the rather slow decay of 
the correlation effect with the distance the polarization of the remaining part of 
the infinite one-dimensional crystal outside the C10H12 cluster amounts to 1% of the 
correlation correction (5IPp°''pj only that is well-below the error bar of the incremental 
scheme itself. 

Let us turn to the electron affinities now. To calculate the correlation corrections 
to the matrix elements corresponding to attached-electron states we proceed as 
described in section 3.4.3. In the case of the diagonal matrix elements EAp°"pj* 
and EAg'^^s*. ,CH2* molecule C10H12 may serve as a minimal-size cluster for the 
correlation calculations. The decay of the second-order increments with the distance 
to the added vr bond-antibond pair of vr symmetry (the most diffuse ones, referred 
to by the vr* orbital only for simplicity) is presented for both matrix elements in 
Table 14.91 The convergence with the distance looks even better than in the case of 
vr hole states, however, that is most probably simply because only a limited number 
of the antibonds are open for excitations from the occupied bonds in the incremental 
scheme for attached electrons which leads to a partial reduction of the virtual space, 
the effect which is more pronounced for increments AEArs{t) with separated bond 
groups r, s and t. 

The convergence of incremental series ()3.34j) with order is shown in Table I^TUI 
the total values of the correlation corrections to the diagonal matrix elements are 
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Table 4.9: Convergence of increments with the distance between correlated bonds 
(in eV). Dashes mean that the increments are too small to be reliable due to size- 
extensivity problem of MRCI. 



^■■^^CH2*,CH2* 


(pi*) 




0.181 


AEA™"pi,(pi*,-l) 


0.314 




(pi*, 


-1) 


0.034 


AEA™rp,,(pi*,-2) 


0.088 


A P A corr 
^^^CH2*,CH2* 


(pi*, 


-2) 




AEA™-,pi,_,(pi*,-2) 





Table 4.10: Correlation corrections to diagonal local matrix elements corresponding 
to antibonds obtained in the cluster C10H12 (in eV). 



matrix element 


AEA=,7() 


E,AEAX(c) 


E,<,AEA:-(c,ci) 




EAcH2*,CH2* 


0.191 


0.816 


-0.161 


0.846 


EApi*^pi. 


0.237 


1.815 


-0.430 


1.623 



presented in the last column of this Table. As in the case of the hole states the 
major contribution to the local matrix elements EA^^ comes from the second-order 
increments AEA^sit) while the third-order contributions AEAj.s{t,u) are at least by 
a factor of 4 smaller than the second-order terms. 

To evaluate the correlation corrections to the off-diagonal matrix elements IPJ;'iT„:, 
and EAq°h2* ch2* 1 '^^ had to use even larger, the C12H14 molecule. The conver- 
gence of the increments with the order and the distance between the correlated 
electrons was checked and turned out to be satisfactory. Summing up these incre- 
ments we get the following values for the corrections: 5EAq°^2* ch2* 1 ~ —0.035 eV 
and 5EA;^°rpi. 1 = 0.365 eV. 

To summarize our results on the correlation corrections to the local matrix ele- 
ments we present in Table I^TTl the SCF matrix elements together with the absolute 
corrections (and their percentage relative the SCF value) and the total correlated 
matrix elements IP'=°'''' and EA'^"" . Note that we also switched back to the notation 
for the matrix elements in the infinite periodic system as done in Tables I3.2H3.5I 
though the correlation contributions are actually calculated in molecules represent- 
ing finite parts of the crystal. In the one-dimensional crystal the polarization outside 
the specified clusters only gives negligibly small effect as compared to the correlation 
effect in the close vicinity of the localized extra charge (see section 4.3) and therefore 
the correlation correction to a local matrix element obtained in the molecules are 
equal to the corresponding one in the crystal within the accuracy of the incremental 
scheme. 

From Table 14.111 one can see that the change of the SCF matrix elements due 
to correlation effects is quite pronounced, ranging from 15% to 30%. One also sees 
that in all cases the absolute values of the SCF matrix elements are reduced, i.e. 
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Table 4.11: Change of the local matrix elements IP i?,nn' and EA/j ^m' of the periodic 
system due to correlation effects. 

matrix element IPp^pipj IPq.chichi IPq.cccc IPp.sigsig 

SCF value 10.526 eV 18.904 eV 21.230 eV 22.396 eV 

reduction -2.007 eV -2.804 eV -3.702 eV -3.867 eV 

19.1 % 14.8 % 17.4 % 17.2 % 

total value 8.519 eV 16.100 eV 17.528 eV 18.529 eV 



matrix element IP±a.pipi IPg.cccHi IPp.cHisig IPq.cc sig 

SCF value 1.683 eV 2.968 eV 3.022 eV 3.022 eV 

reduction -0.316 eV -0.718 eV -1.040 eV -0.936 eV 

18.7 % 24.2 % 34.4 % 31.0 % 

total value 1.367 eV 2.250 eV 1.982 eV 2.087 eV 



matrix element EAp^pj* pj* EAq^qhi* chi* EA^a^pj* pi» EA^a.cHi* chi* 

SCF value -4.497 eV -5.603 eV -1.054 eV 1.012 eV 

reduction 1.623 eV 0.846 eV 0.365 eV -0.035 eV 

36.1% 15.1 % 34.6 % 3.5 % 

total value -2.941 eV -4.764 eV -0.689 eV 0.976 eV 
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overall the ionization potentials decrease while the electron affinities increase. This 
indicates that the reduction of the total energy due to electron correlation is always 
stronger in the charged systems than in the neutral one. 

4.2 Correlated band structure of ^rans-polyacetylene 

Having obtained the local matrix elements including all correlation effects (the 
"total values" in Table I^TTjl we are able to compile the band structure of trans- 
polyacetylene which takes into account electron correlation, the ultimate goal of our 
approach. For this purpose we replace the LMEs in the truncated series ()2.32|) and 
fl2.33|) obtained on the SCF level by the correlated ones and diagonalize the obtained 
matrices at each /c-point analogously to the case of the SCF bands. The new band 
structure (we call it correlated band structure) will exhibit small deviations from 
the one obtained by an infinite summation in Eqs. ()2.32|) and ()2.33p like the SCF 
band structure (see Figs. ITBl and IXT^ . To correct for this (minor) deficiency we 
simply add at each /c-point the difference between the band energies obtained by 
the local SCF matrix elements and those obtained by the CRYSTAL code for an 
infinite periodic tPA single chain. The resulting correlated band structure of trans- 
polyacetylene is presented in Fig. 14.31 (red lines) and compared to the SCF band 
structure by the CRYSTAL code (black lines). 

From this Figure we see that the band structure changes dramatically when 
electron correlation is taken into account. Two effects are clearly seen: i) a shift of 
the " center-of-mass" positions of all bands (valence bands are shifted upwards and 
conduction bands downwards) and ii) a flattening of all bands (signiflcantly more 
pronounced for the valence bands). These effects lead to a narrowing of the band 
gap, a change of the sign of the electron affinity of trans-polyacetylene and an overall 
reduction of the band widths. But there are also some more subtle changes like the 
narrowing of the avoided crossing of the 2nd and 3rd lowest valence bands close to 
the X point or the lifting of the accidental degeneracy of the highest two valence 
bands at the F point. 

The value of the band gap of the infinite tPA single chain is reduced by 2.31 
eV from 6.42 eV (on the SCF level) down to 4.11 eV when electron correlation is 
included, a 36.0% reduction. Similar values of the reduction of the SCF band gap, 
2.35 and 2.38, are reported in two recent independent MP2 studies, |SB96j and 
AKSOlj . respectively, on the tt bands of isolated tPA chains. 

The corrected values of the lowest ionization potential and the highest electron 
affinity of trans-polyacetylene are 4.58 eV and 0.47 eV, respectively. The posi- 
tive sign of the electron affinity (after inclusion of electron correlation) tells that 
trans-polyacetylene is able to retain an extra electron as is known from experiments 
IKP+891 . 

The changes of the band widths are summarized in Table 14.121 There, the 
a valence band width refer to the width of all four a valence bands while the a 
conduction band width is that of the two a conduction bands only. One sees that 
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Figure 4.3: Correlated band structure of the infinite trans-polyacetylene single chain 
(red lines) compared to the SCF one (black lines). 

Table 4.12: Change of the band widths due to the correlation effects (in eV). 



reduction 



bands 


SCF width 


correlated width 


absolute 


relative 


a valence 


17.30 


13.01 


4.29 


24.8 % 


TT valence 


7.40 


6.14 


1.26 


17.0 % 


TT conduction 


4.86 


4.67 


0.19 


3.9 % 


a conduction 


5.25 


5.05 


0.20 


3.9 % 
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the relative reduction of the conduction bands widths is much smaller than that of 
the valence bands. Since the dispersion of the bands is mainly determined by the 
hopping matrix element between nearest-neighbor sites (especially in ID systems) 
one can conclude that the hopping of holes is much stronger affected by electron 
correlation than hopping of extra electrons. 

4.3 Accuracy of the method 

Now that we have obtained numerical results for some physical quantities which 
are characteristic for the investigated system the question arises how accurate these 
results are and whether we are able to estimate the accuracy of our method. In 
this section we show that we can indeed provide some estimated error bars for the 
calculated quantities. 

As the correlation corrections to the local matrix elements (the key quantities 
of our method) are composed of dozens of increments we first estimate the error 
for the individual increments. As can be seen from equations fl3.15|) and ()3.2U|) - 
fl3.24|) each individual increment is defined as the difference of several correlation 
energies (correlation corrections to the total energies of various states under various 
conditions which are calculated to evaluate a particular increment). The uncertainty 
in the evaluation of these individual increments essentially arises from the size- 
extensivity error of the CI(SD) method. As shown in section 3.5 we are able to 
reduce this error to about 10~^ of the value of the correlation energy (see last 
columns in Tables ITHl and ITTj) . In all our calculations these correlation energies are 
of the order of a few electron-volts. Therefore, we estimate the error for a single 
increment being a few meV. That is why in our calculations we simply omit all 
increments whose absolute values are smaller than 1 meV (see Tables I^TI and . 

The possible accumulated size-extensivity errors in the correlation correction to 
some particular local matrix element are expected to be not higher than the error 
for a single increment times the number of increments entering the series. Usually 
this amounts to few dozens of meV. 

Another kind of error in the correlation corrections to LMEs emerges from the 
truncation of the incremental series after some order (after the third order in our 
case). In section 4.1 we have explored the accuracy of the truncated incremental 
scheme in detail. In the case of the most compact localized orbitals {a bonds) the 
error produced by the incremental scheme is of the same order as the accumulated 
size-extensivity error (4-7% of the total value of correlation correction to a LME, 
see the discussion of Tables Wl] and . However, in the case of the more diffuse 
TT orbitals the error due to the truncation of the series (~15%, see the discussion 
of Table 14. 5j) is higher than the estimated size-extensivity error for diagonal matrix 
element (~2%). This indicates that for more diffuse orbitals the fourth term of 
the incremental series ()3.18p may give noticeable contribution. Unfortunately, we 
can not calculate this term explicitly as individual increments of the fourth order 
are usually smaller than their size-extensivity error. Thus, we agree that the er- 
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TOT for correlation corrections to the diagonal matrix elements IPpi,pi, EApj* pi* and 
EAcHi*,CHi* are not more than 15% and those for all the other matrix elements not 
more than 10% of absolute values of the corrections. 

Another question which we would like to discuss here concerns the effect of po- 
larization of infinite insulator due to the presence of an uncompensated charge. The 
polarization of the infinite crystal outside the explored cluster may also contribute 
to correlation corrections. We can estimate such a contribution by a continuum 
approximation as already done earlier by other authors, e.g. in |HHF84j . |GSF97j 
and [AFSnnj . 

In this approximation far-range contributions to the reduction of the SCF energy 
of the (A^ — 1)- or (A^ + l)-electron system are attributed to the polarization energy 
of an insulator in an electric field. This energy is equal to 

SEf"' = ^JpEdV (4.1) 

where E is the electric field, P = -^^^E is the polarization response and e is the 
dielectric constant of the crystal. The subscript i in Eq. ()4.1|) denotes that this 
energy corresponds to the charged system with a hole (an extra electron) being in 
the bond (antibond) i. In other words, we are estimating the polarization effects on 
diagonal matrix elements here. 

Let us say that starting with a distance ri from the center of a localized orbital 
we can treat the electron in this orbital as a point charge. Then the integration in 
Eq. ()4.H) can be divided in two parts: 

SEr' = 5i5.(ri) + ^"dr (-^)' = 6E,{n) + CrY\ (4.2) 

Here we performed the integration for the case of one- dimensional system, 6Ei[ri) 
is the change of the energy due to electron correlation for distances less than ri from 
the localized charge (this value we calculate explicitly by the incremental scheme) 
and C is a constant characteristic for the given material. Knowing two SEi{r) values 
SEi{ri) and SEi{r2) for two different cut-off radii ri and r2 (r2 > ri) we can extract 
the tail correction Cr^^ and by this estimate change of the energy of the whole 
charged system 

e p . ^ 6E,ir,)ir,/r,r-6E,in) 

^E.ioo) = ^-7^7)331 • (4.3) 

Note the different power (r^) which enter this relation compared to the corresponding 
relation for a bulk crystal |AFSnnj r^. 

Thus, the correlation correction to diagonal matrix elements with the polariza- 
tion of the infinite crystal taken into account is 

ATP— _ SiP:Tir2){r,/nf-5iP:Tin) 

^ '''' ~ {r^/nf - 1 ^ ' 
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where (5IP^°"(rj), i=l, 2 are the correlation correction to the matrix element IPaa 
when all bonds up to distances rj from the bond a are correlated. An analogous 
formula holds for 5EA^°"(cx)). 

The hopping matrix element IPq;, is one-half of the difference between energies 
of the states ($a + ^i,)/\/2 and ($a — $6)/-\/2 which are two states with holes of 
two different shapes delocalized within two bonds a and b. Therefore, the change 
of the hopping matrix element in continuum approximation emerges due to the 
polarization of the insulator in the electric field of a dipole. Substitution of the 
field of a point charge in Eq. (j4.2|) by that of a dipole leads to the approximate 
formula for the correlation correction to the off-diagonal matrix elements with the 
polarization of the entire crystal taken into account: 

^IP^^ (^) = (r,/r,)^ - 1 • ^'-'^ 

Let us calculate explicitly the missing polarization effect in the correlation cor- 
rection (5IPpj pi caused by adding farther n bonds using the data on the decay of the 
second-order increments in Table 14.71 The total value of the correlation correction 
in a cluster C18H20 (four unit cells on both sides from the central one) as obtained by 
the incremental scheme effectively is 2.007 eV and the cut-off radius r2 = 4a. As the 
correlation effect associated with the a bonds decay much faster than those arising 
from the vr bonds, the change of the correlation correction when the effective cluster 
increases from C10H12 to C18H20 is mainly due to correlation of the electrons in it 
bonds. Thus, the value 51?™ pi(2a) is equal to -1.871 eV. Using Eq. (jO)) we obtain 
the value of the correlation correction in the infinite system (5IPpj'j^i(oo) = 2.026 eV 
compared to 51?™ pi(4a) = 2.007 eV. Therefore, the polarization of the outer part 
of the crystal leads to less than 1% change of the correlation correction which can 
be neglected within the error bar of the incremental scheme. 

The continuum correction in the case of the off-diagonal matrix elements is 
expected to be even less as the decay of correlation effect with the distance is more 
rapid than for diagonal ones (see Eqs. 14.31 and 14. 4j) . In fact for an effective CigHig 
cluster (with r = 3.5 a) 5IPp°pj ^^(3.5 a) = —0.316 eV while for the smaller effective 
CgHio cluster (with r = 1.5a) 51?™ pi^^(1.5a) = —0.148 eV. The resulting continuum 
correction is -0.002 eV leading to 51?™ pj^(oo) = —0.318 eV. Thus, we conclude that 
the effect of polarization of the infinite one-dimensional crystal outside the treated 
cluster is negligibly small. 

Now we can estimate the accuracy of our method for the ultimate physical quanti- 
ties characterizing the investigated material. The shift of the ionization potential due 
to electron correlation is governed by two corrections: = 5IPQ°p[pj — 251P™p[pj. 

This is easy to realize in terms of a Hiickel (or tight-binding) model where only 
diagonal matrix element and the hopping between nearest-neighbor sites are con- 
sidered. Then, energy of the valence vr band as a function of crystal-momentum 
k is e^{k) = -IPo°p[,pi - 2IP'i°^[pi cos(aA;). The top of this band is at the X point 
{k = 71 /a) and therefore the ionization potential is given by IP = IPo°p[pj — 2IP™p[pj. 
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We have estimated the error for the correction 5IPQ°p[pj being equal to 15% of its 
value (-2.0 eV) which is 0.3 eV. The error for the correction 5IPQ°p[pj is smaller by 
one order of magnitude and is omitted. Thus, we consider the ionization potential 
of tPA single chain to be IP = 4.6 ± 0.3 eV. Approximately the same error bar is 
obtained for the electron affinity: EA = 0.5 ±0.3 eV. Consequently, the error for the 
band gap is at most 0.6 eV that is 15% of its value. However, we have to emphasize 
that the estimated error of 0.3 eV is an upper bound for the energy of any band at 
any fc-point which is in the order of 1% of the total range of band energies only. 

As one can conclude, the two main reasons, which lead to the error bar for 
correlation corrections up to 15% of their values, are i) the error due to truncation 
of the incremental scheme after the third order and ii) the size-extensivity error of 
the CI(SD) method for evaluation of correlation energies of (A^± 1)- and A^-electron 
states. The truncation error of the incremental scheme can not be reduced by taking 
into account the next order because the values of the increments of the fourth order 
are below the estimated size-extensivity error for the entering correlation energies. 
Therefore, namely the latter error controls the accuracy of our whole approach. 
To improve the accuracy further one first needs a true size-extensive multireference 
correlation method that could be multireference coupled-cluster which is at the stage 
of development at the present moment. 

4.4 Comparison to experimental data 

In this section we would like to compare our result for the band gap to experimental 
data. One has to emphasize that there are, of course, no experimental data on 
trans-polyacetylene single chains and all data available correspond to bulk trans- 
polyacetylene and we have to compare properties of a bulk system with calculated 
results for the linear one. 

One attempt to estimate the effect of the interchain coupling on the band gap 
of tPA single chains was made by Ayala et al. |AKS01j . They have compared the 
band gap of an isolated single chain with the band gap of a system consisting of two 

o 

chains one precisely on the top of another at a distance of 3.335 A. In this particular 
arrangement of the two chains the intermolecular overlap of the vr orbitals (which 
define the gap) is maximal. The chosen distance corresponds to the minimum of 
the interaction energy of the two stacked chains on the MP2 level of approximation. 
The reduction of the band gap due to interchain coupling was found to be 1.22 eV 
on the HF level and further a 0.22 eV reduction (1/6 of the HF value) was obtained 
when electron correlation was taken into account (the M0ller-Plesset perturbation 
method to second order was used for this purpose). Thus, one can conclude that the 
main reduction of the band gap due to the interchain interaction already emerges 
on the HF level. 

We investigate this reduction for the true fishbone-like packing of tPA chains in 
the bulk material as shown in Fig. 13.21 which was found experimentally |FC+82| . 
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The unit cell of bulk tPA consists of two C2H2 units of linear tPA (one on each non- 
equivalent chain). The interchain interaction leads to the splitting of the otherwise 
two-fold degenerate bands of single tPA chain. This splitting at the X point amounts 
to a 0.6 eV reduction of the band gap of the bulk system as compared to that of a 
single chain. Further reduction of the band gap due to electron correlation in each 
of the four nearest-neighbor chains is not expected to be noticeable because in the 
experimentally-found packing of the chains the overlap of the it orbitals between 
neighbor chains is actually minimal. Thus, based on the MP2 data presented in 
|AKS01j we estimate it to be less than 1/6 of the HF reduction for each of the four 
chains leading at most to a further 0.2 eV narrowing of the band gap. Therefore, 
the account for the interchain splitting reduces the calculated band gap of the bulk 
tPA down to 3.3 eV. 

Two other effects should also be mentioned which can lead to a reduction of the 
calculated band gap of tPA. Firstly, though our basis set (VTZ) is already rather 
flexible we have not yet reached the basis set limit. A more rich basis set will 
affect the correlation corrections of the LMEs as a larger variational space will be 
provided to the MRCI calculations. Secondly, the long-ranged polarization of the 
infinite bulk system has to be considered. There, one expects a more pronounced 
effect than in the case of a one- dimensional system. A discussion of the polarization 
effect on local matrix elements of different bulk systems can be found in |GSF97j 
and |AFS00j . However, both effects mentioned here are small compared to the 
calculated electron correlation effects and each of them can lead to a few tenth of 
an electron-volt reduction of the band gap only. 

Next, we would like to address the experimental value of 2 eV for the tPA 
band gap obtained from absorption spectra of bulk trans-polyacetylene | FO-|-79 |, 
TG-l-801 and |L88j . There is no clear indication in these publications that the 



observed maximum of the absorbtion coefficient at approximately 2 eV can really be 
assigned to the interband transition rather than to an excitonic state. A theoretical 
investigation of such gap states showed |RL99j that the energy of these excitons in 
tPA show up approximately 0.4 eV below the band gap. Therefore, the fundamental 
band gap of bulk trans-polyacetylene is expected to be 2.4 eV. 

Thus, taking into account the bulk effects in tPA which all lead to a reduction 
of the band gap with respect to the one in a ID single chain and also the estimated 
error of our method the calculated gap reaches the experimental value. 
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Chapter 5 
Excitons 



In this Chapter we would hke to show how our approach to correlated band struc- 
tures can be extended to the study of the electron correlation effects of neutral exci- 
tations in periodic systems, namely excitons, such that standard quantum-chemical 
methods can be used. We refer here to the process of the absorption of a light 
quantum by the electrons of an extended non-conducting system in its ground state 
i.e. the process of generating an excited state of the same neutral system. In this 
case no electron leaves the system (in contrast to the photoionization process which 
we refer to when describing the valence bands). 

In frames of the most simple model one can say that an electron from the valence 
band u and with band energy e^i, is excited to the conduction band fi where it 
occupies the virtual Bloch orbital ipk'fi with energy Ek'n- The wavefunction of such 
a state can be obtained by applying the two operators Cku and cj,,^ entering the 
equations ()2.22|) and ()2.23|) to the ground-state wavefunction ()2.21|) of the closed- 
shell system (we omit here the spin indices as we suppose that the electron does 
not change its spin during such a process). The operator Cku produces a hole in the 
valence band ly and the operator c\,^ creates an electron in the conduction band /i. 
If we neglect the electron-hole interaction, the energy to produce shuch an excited 
state is equal to Ek'^ — £ku- This way we describe a pure interband transition and 
the threshold energy for such processes is equal to the fundamental gap of the non- 
conducting system. 

However, if we take into account the interaction between the electron in the 
conduction band and the hole in the valence band, the system may gain energy 
from the electron-hole interaction when a bound electron-hole pair is formed. In 
other words, in the extended system there exist gap states which correspond to ex- 
citons (bound electron-hole pairs) travelling through the crystal. Such gap states 
were found experimentally in many semiconductors and insulators, in particular in 
trans-polyacetylene |LE+81| , |OB82j and |SY-f-82| , the system which we refer to 
here. These states form bands and we would like to use our approach to calculate 
these excitonic bands. This means, we want to construct a set of localized model 
electron-hole pair wavefunctions (forming our excitonic model space) which can be 
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represented in finite clusters, and calculate the corresponding local matrix elements 
(diagonal and off-diagonal) in sufficiently large clusters of the crystal. These exci- 
tonic LMEs are essentially the matrix elements of an effective Hamiltonian obtained 
from an equation similar to Eq. ()2.(j2j) . and we will get the excitonic bands by the 
use of a (back) transformation analogous to Eqs. ()2.32|) and ()2.33|) . 

We exclusively refer here to bound excitons which carry a single crystal momen- 
tum (associated with the center of mass motion) and are otherwise regarded as a 
fixed entity. The appealing aspect of this approach is that in quantum chemistry 
which deals with finite systems an exciton is just an eigenstate of the given system 
like the ground state, but different from it, and the same CI ansatz and methodology 
as used for the correlated ground-state wavefunction can be used. In fact, the CI 
equations are the same and one is just looking at another root of these equations. 

To describe a bound electron-hole pair of a crystal in real space it is convenient 
to use local Wannier orbitals (both occupied and virtual ones) to construct the 
local excited configurations. Let as consider one such configuration where an 

electron is removed from the occupied Wannier orbital ipB.n centered at the unit 
cell with lattice vector R and put into the virtual orbital (pn'm centered at the unit 
cell with lattice vector R' while the rest of the system remains unchanged (i.e. we 
start again from the frozen orbital approximation). Taking into account Eqs. (j2.15|) - 
fl2.19|) one can readily show that the energy to produce an excited state with this 
configuration is equal to 

^exc ^ ^^R^m I ^ I ^R^m^ _ 

= Em — £n — {{R'm Rn I R'm Hn) — (R'm Rn | Rn R'm)) 

= -EA^ + IP„- ((R'mRn|R'mRn) - (R'mRn|RnR'm)) (5.1) 

where we have used | Rn) instead of | (p-Rn), and Em and En are the energies of 
the Wannier orbitals (associated with the bands fi and u, respectively) which are 
independent of the vectors of lattice translation R' and R. The energy E^^^ is 
not equal to the difference of the ionization potential IP„ and the electron affin- 
ity EAm which defines the interband transition but differs from it by the term 
(R'm Rn | R'm Rn) — (R'm Rn | Rn R'm). This term decays with increasing dis- 
tance between the centers of the Wannier orbitals ipn^ and v^R'm (as follows from 
Eq. (|2.6|) ) and is the interaction energy between the electron and the hole. Such 
terms will be responsible for the formation of bound electron-hole pairs when we 
describe excited states of crystals (or big clusters) as linear combinations of local 
configurations like ^rjJ". The energies of such states can be substantially reduced 
compared to the energy of the interband transitions and thus the gap states of 
semiconductors and insulators may emerge. 

Below, for the sake of simplicity, we consider a single valence and conduction 
band only that yields one occupied and one unoccupied Wannie orbital per unit 
cell. Such a simplification is justified for irans-polyacetylene, when excitations of the 
valence tt electrons to the lowest conduction band of tt symmetry are considered. We 
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then can omit the band indices n and m baring in mind that subscripts (superscripts) 
in correspond to the hole (electron) orbitals. Also, again we will refer to a 

one-dimensional system as the excitons in the conjugated vr system of polymers are 
essentially confined to a single polymer chain |FG-|-99| . 

Let us start to describe a bound electron-hole pair by writing an excitonic wave- 
function in the form of a linear combination of local configurations $^ : 

^^^^ = E ^rW^R = E E ^M^l^' (5.2) 

RR' R d 

where we introduced the relative distance d of the centers of the Wannier orbitals 
and ipRi: d = R' — R. Here we have to make a few remarks. Let us focus for a 
moment on a wavefunction where the hole is pinned to a certain unit cell, i?, say. 
Then the CI ansatz ()5.2|) for the exciton is reduced to 

d 

Firstly, if we set all coefficients aR^d ^ 0) to zero and ^^^(O) = 1 we get a so-called 
Frenkel exciton |Fr31j where the excitation is restricted to one unit cell. Yet, we use 
the ansatz which allows us to describe excitons in the general case when the 
atomic structure of the given crystal is explicitly taken into account. Secondly, it is 
possible to produce A^^'''* excitonic states using the ansatz ()5.3|) where N™^ is the 
number of virtual WOs in the system. However, we focus our attention on the ener- 
getically lowest excited state which corresponds to the gap state and is usually well 
separated from the other excitonic states by some energy. Thirdly, doubly excited 
configurations like <J)^^^'^'^+'^ may, of course, also contribute to the wavefunction 
of the excited state but here we aim to construct our model space, which is built 
up by singly excited configurations which describe the electron-hole pair as being a 
two-particle object. The effect of doubly (and higher) excited configurations on the 
energy of excited states can be taken into account later by the MRCI procedure. 

The representation of the exciton-state wavefunction given by the right-hand- 
side expression in Eq. (j5.3|) suggests to plot the distribution of the electrons relative 
to the fixed hole. As the electron-hole interaction energy (the two last terms in 
Eq. (jOJ)) contained in each configuration decays with the distance between the 
electron and the hole orbital, one expects configurations with small distances d to 
dominate the wavefunction ()5.3|) . This expectation is based on the fact that the 
electron tries to stay as close as possible to the hole. 

An investigation of such electron-hole pairs in a trans-polyacetylene single chain 
was presented in Ref. |RL99j which was done in the framework of a different (re- 
cently developed) approach based on solving the Bethe-Salpeter equation of the 
two-particle Green's function and using the GW approximation for the self-energy 
operator. That approach allowed the authors to visualize the distribution of the 
electron relative to the hole fixed somewhere in real space (see Fig. 15. ip . From this 
Figure one concludes that such localized objects can indeed be found in periodic 
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Figure 5.1: Electron-hole wavefunction of the exciton state of a tPA single chain 
showing the distribution of the electron relative to the hole which is fixed at some 
arbitrary position on the chain (at A). The figure is taken from Ref. jRL99j where 
the calculations were done within the GW approximation and using the Bethe- 
Salpeter equation. 



non-conducting systems and the probability to find the electron of this pair at some 
place decays with the distance from the fixed hole. Hence, a finite cluster can be 
used to represent this localized excitonic state. 

As follows from Fig. 15.11 the dimension of a single electron-hole pair in a tPA 
chain is approximately 30 unit cells that implies that a CgoHgo cluster is required 
to properly describe such an exciton within our cluster approach. Presently we can 
not stand calculations for such big clusters of irans-polyacetylene with reasonable 
basis sets. Therefore, all our calculations presented in this Chapter are done for a 
much simpler system consisting of 9 H2 molecules which are arranged analogously 
to the system of 4 H2 molecules considered in section 3.5 (see Fig. I3.14ll but with the 
distance of 3 A between neighbor H2 molecules to allow the exciton to travel along 
the chain. This system represents a 9-unit-cell cluster of an infinite (H2)a; chain 
which can be considered as a simple model for the vr system of a tPA chain. In our 
9 X H2 cluster we first generate localized occupied and virtual molecular orbitals (in 
fact, bonding and anti-bonding orbitals of the H2 molecules) and then compute the 
energetically lowest excited state when only the configurations $5^*^ with the hole 
fixed at the central molecule H2 (the 5-th one) and the electron in any of anti-bonding 
orbitals in the 9 H2 chain are used (by a reference space only Cl-calculation using 
the MOLPRO program package |M2nnnj . |WK88j . |KW88j . |KW92j ). The squares 
of the coefficients a^ld) of these configurations, being the probabilities to find an 
electron in the localized virtual orbital centered in the unit cell R' = 5 + d, are 
presented in Table 15.11 We find the same type of electron distribution around the 
fixed hole as reported for a tPA single chain in Ref. |RL99j (see Fig. l5.ip . Therefore, 
the electron-hole interaction, which enters the state energy via terms analogous to 
the two last terms on the right-hand side of Eq. 15. H leads to the formation of well- 
bound electron-hole pair and such a (pinned) pair in a crystal can be described in 



97 



Table 5.1: Sqares of coefficients 05(6?) of the configurations ^^'^'^ entering the local 
excitonic wavefunction (j5.3|) in a 9 x H2 cluster of an (H2)a: chain with the hole being 
fixed at the central (5-th) H2 molecule. 



d 





±a 


±2a 


±3a 


±4a 


(«5(rf))^ 


0.951 


0.024 


1.54 X 10-^ 


5.84 X 10-^ 


1.79 X 10-9 



the form of the linear combination ()5.3|1 with only a finite amount of local excited 
configurations 

Such localized states given by Eq. (j5.3p are not eigenstates of the periodic system 
as they are degenerated. To remove the degeneracy we proceed as in chapter II of 
Ref. |K63j that is the pinned-hole excitonic wavefunctions are treated as rigit 
objects and a set of LCAO ansatz is made for the wavefunction of the excited state 
in the one-dimensional crystal (like in a tight-binding model) : 

R R d 

Here K is the crystal momentum of pair made up by the hole and the electron being 
distributed among virtual WOs in the close vicinity of this hole. The ansatz (j5.4j) 
takes into account the periodic structure of the crystal and the excited states of 
the crystal can be referred to "particles" of neutral electronic excitations (which are 
called excitons) travelling through the crystal with wave vector K. The energies 
to produce such excited states regarded as a function of K form the excitonic 
bands which can be obtained from the finite-cluster calculations. 

In the 9 X H2 cluster we open all 9 bonding and 9 anti-bonding orbitals to 
construct the singly excited configurations $^ which form the model space and 81 
excited states can be obtained as linear combinations of these configurations: 

9 9 

a=l s=l 

where index a (s) denotes an occupied (virtual) localized molecular orbital. Actu- 
ally this equation is just the finite-cluster version of the crystalline ansatz given in 
Eq. ()5.2|) . Among all these excited states we find 9 states which are energetically 
well-separated from the rest and lower than all other excited states. These 9 states 
i = 1, . . . , 9 are strongly-bound excitons we are looking for here and can be 
regarded as some linear combinations of 9 suitable localized excited states. These 
9 "delocalized" states E^^'^ are easily identifiable since the energies of these states 
only spread over an energy interval of 1.4 eV (this value is approximately equal to 
the excitonic band width) while they are separated by a gap of 2.1 eV from all the 
other excited states. In fact, we focus our attention here on the lowest excitonic 
band which corresponds to the "ground state" of an exciton. 
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To calculate the excitonic band in the framework of our approach we have to 
define local excitonic matrix elements (which are stransferable from clusters to the 
infinite periodic system): 

EX,, = - 6,,Eo = ($r I H I $f ^) - 6,,Eo (5.6) 

in terms of some /oca/uec? excitonic wavefunctions ^f^^. The latter are obtained by a 
unitary transformation applied to the A^"'^'^ delocalized excitonic wavefunctions "^f^^ 
in the cluster 

where the unitary matrix uji is determined by minimizing the total spatial spread 
of the electron-hole pairs 

1=1 

(cf. Eq. ()3.5|) ). Below we provide one possible definition of the spatial spread of the 
excitonic wavefunction in terms of singly excited configurations. 

Note that the introduction of an electron-hole pair with a pinned hole was only 
used to develop the idea of the quantum-chemical description of a well-bound exciton 
in terms of singly excited configurations. However, in both, crystals and clusters 
there exist no such fixed hole, and electrons and holes together are delocalized. 
Hence, it seems natural to describe the bound electron-hole pairs in terms of exci- 
tonic coordinates, i.e. a coordinate R of the "center of mass" of the exciton and a 
relative coordinate f. They can be defined as the population-weghted average of the 
"centers of mass" and the relative coordinate R^ and of each of the configurations 

entering the excitonic wavefunction: 

a,s a,s 

(cf. Eq. dSHI)) where 

K = liRa + Rs) and r', = ^{Ra-Rs), (5.10) 

Ra denotes the center of the bonding orbital and Rs the center of the anti-bonding 
orbital ips- 

The spread of each electron-hole pair is defined by the average distance between 
the electron and the hole (given by f,) and by the average deviation of the "center 
of mass" of each contributing configuration R^ from the " center of mass" Ri of the 
exciton itself: 



Y.{K-R^f{a,lY (5.11) 
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(cf. Eq. (jSISD). 

Having obtained the unitary matrix Uji which corresponds to the minimum of 
the functional P^"^ in Eq. (|5.8p or some other measure for the extent of the exciton, 
we get A^*^^^ localized excitonic wavefunctions Then, the matrix elements of 

the Hamiltonian H^^'^ can be calculated analogously to the elements of the effective 
Hamiltonian for hole or electron attachment states as given by Eq. ()2.62|) : 

i' 

In the periodic system one expects to find one such localized excitonic wave- 
function ^^'^ in each unit cell (with R being the translation vector of this unit cell) 
which all together correspond to the lowest excitonic band of the crystal. Note, that 
here the center of the exciton is its " center of mass" Rr in contrast to the position 
of the fixed hole in the ansatz ()5.3p . The local excitonic matrix elements in the 
periodic system are then defined as 

Ex^ = - ^orEo = ($r I H I $r) - ^o^^^o (5.13) 

The diagonal matrix element EXq is the (Brillouin zone averaged) energy of the 
exciton and the off-diagonal elements EX/j are the hopping matrix elements between 
the reference (zero) unit cell and the unit cell at R. As the cluster represents a 
finite part of the crystal one can associate the excitonic LMEs EX,y obtained in the 
cluster (for the two localized excited configurations and •l'^™ with centers Ri 
and Rj, respectively) with the crystalline LMEs EXj^_^^. The closer the centers of 
the chosen localized excitonic wavefunctions ^i^'^ and ^T'^ are to the center of the 
cluster the better the LMEs EXj^ from the cluster approximate the corresponding 
crystalline LMEs EXr. 

The excitonic band finally can be obtained by the local excitonic matrix elements 
(j5.13|) being subject to a transformation similar to Eqs. (j2.32|) and (|2.33|) : 

EX(A;) = ^ e^'^^EXjj. (5.14) 

R 

As the local off-diagonal matrix elements EX/j are expected to decay with increasing 
distance R between the centers of the localized electron-hole pairs <I>q™ and the 
series (I5.14p can be truncated at some cut-off radius -Rcut such that the omitted part 
of the series is negligible with respect to the terms already summed up. Therefore, 
one only needs a finite number of the off-diagonal matrix elements to calculate the 
excitonic band with a given accuracy, and both, diagonal and off-diagonal excitonic 
matrix elements can be obtained from finite clusters. 
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Chapter 6 

Conclusions and perspectives 



The main result of the present thesis is the implementation of a multireference 
configuration interaction method with singly and doubly excited configurations to 
account for the electron correlation effects in both neutral and charged infinite peri- 
odic systems, in particular in excited electron hole and attached electron states. Our 
approach is one of the first approaches for solids which is purely based on the cor- 
related wavefunction of the system. Sufficiently large basis sets (valence triple zeta) 
to describe the correlations especially in the anionic (A^ + l)-electron systems are 
used. As an output we obtain the quasiparticle band structure (all valence and the 
lowest conduction bands) of a non-conducting material where electron correlation is 
systematically taken into account which allows to extract the ionization potential, 
the electron affinity, the band gap and the band widths of the material as those 
properties which can be measured in experiments. All these quantities are obtained 
with controlled accuracy. Thus, the approach provides highly desired information 
on electronic properties of extended systems on a very high level of sophistication 
and, as we believe, will find many interesting applications. We have applied it to 
irans-polyacetylene single chains and obtained complete and quantitatively correct 
information on electronic properties of the investigated system. Further, an ex- 
tension of the above approach to excitons (i.e. optical excitations) in crystals is 
developed which allows to use standard quantum-chemical methods to describe the 
electron-hole pairs and to finally obtain excitonic bands. 

Our method starts from the many-body wavefunction of the infinite closed-shell 
system obtained in the Hartree-Fock approximation. Local quantities in real space 
(namely, localized electron orbitals and local matrix elements) are then defined 
which are well-transferable from finite clusters to the infinite system. By the use 
of these local matrix elements derived from clusters the SCF band structure of the 
extended system is accurately reproduced. Exploiting the local character of the 
electron correlation we perform correlation calculations by the MRCI(SD) method 
in the finite clusters for the A^-, (A^ — 1)- and (A^ + l)-electron states and as a result 
obtain the change of the local SCF matrix elements due to the correlation effects 
which again turned out to be well-transferrable from the clusters to the periodic 
system. 
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To reach the current stage of implementation of our approach the following 
challenging problems were solved. A proper treatment of the (A^ + l)-particle states 
requires localized virtual SCF orbitals which are essentially the same for both infinite 
periodic systems and clusters. Such orbitals can not be generated in a cluster by 
any standard localization scheme applied to some ad hoc set of the canonical virtual 
orbitals of the cluster (in contrast to the occupied orbitals). Therefore, we first 
generate a set of localized virtual Wannier orbitals in the crystal which correspond 
to bands of our interest and then project them onto the space of the virtual canonical 
orbitals of a proper cluster. This way we are able to generate localized virtual 
molecular orbitals and the relevant local SCF matrix elements in the clusters which 
are used further in the correlation calculations. 

To systematically account for the correlation effects the incremental scheme was 
exploited. However, this scheme (originally designed for the ground and hole states 
calculations) has the problem of the emergence of a divergently large number of 
increments when bigger clusters are considered. To handle this problem the scheme 
for the hole states was reformulated in terms of orbital groups. Also, the scheme 
was extended to the case of attached-electron states. 

We showed that within the multireference configuration interaction method one 
can also handle the problem of satellite states. This problem appears in systems 
with a relatively small band gap and leads to severe instabilities in the calculations 
when many-body perturbation theory is used to account for the correlation effects. 
In the MRCI method the satellite states can be treated together with the true hole 
or attached-electron states by a proper modification of the reference space however 
by the cost of increasing computational effort. 

The use of configuration interaction methods with singly and doubly excited 
configurations in combination with the incremental scheme requires a very accurate 
handling of the size-extensivity error intrinsic to all truncated CI methods. We 
could show that the known Pople correction applied to the correlation energy of 
the neutral closed-shell system reduces this error to the required level. However, 
there did not exist any correction for the size-extensivity error of the MRCI(SD) 
method, which in case of the charged open-shell systems could lead to a comparable 
precision. In the present work we have developed such a correction. The analytic 
formula for it is derived and its performance is checked. The new correction allowed 
us to extract the final results with a reasonable accuracy. 

The whole approach is focused on obtaining local matrix elements which cor- 
respond to correlated wavefunctions of the system and contain information on the 
N-, [N — 1)- and {N + l)-electron states of a system when electron correlation is 
systematically included. Using these matrix elements we were able to compile a 
band structure (both valence and conduction bands) of the chosen crystal which 
fully takes into account correlation effects. From this band structure we extract 
correlation-corrected values of some characteristic electronic properties of the sys- 
tem which can be measured experimentally. They are the band gap, the band 
widths, and the ionization potential and electron affinity. We applied our method 
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to an infinite irans-polyacetylene single chain. We found that upon inclusion of 
electron correlation, the calculated results approved substantially as compared to 
those on the HF level and that they are in good accordance with the experimental 
data once the effect of interchain interaction and long-range polarization in bulk tPA 
are also taken into account. 

As a remarkable feature of our method we emphasize the possibility to estimate 
the error for the obtained quantities and to systematically improve the results to 
any desired accuracy by just increasing the computational effort. For the calculated 
correlated band structure of a irans-polyacetylene single chain which spans an energy 
range of approximately 30 eV the estimated error for any band energy at any /c-point 
does not exceed 0.3 eV. 

Comparison of our results on a tPA single chain to the experimental data on 
the bulk crystal, composed of weakly interacting chains, indicates that the study of 
correlation effects in polymers can not be reduced completely to the treatment of 
isolated chains. There is clear evidence that accounting for the mutual interaction 
of the tPA chains in the bulk material improves the theoretical results with respect 
to the experimental data. The correlation of electrons in neighbor chains has to 
be considered for that purpose. This can be done with our approach. However, to 
calculate these effects explicitly is beyond the scope of the present study. 

Finally, a new wavefunction-based approach which allows to describe well-bound 
excitons in crystals by means of standard quantum-chemical methods is established. 
Exploiting the finite dimension of a bound electron-hole pair (extending over several 
unit cells) we can restrict ourselves to finite clusters of a crystal to calculate the 
relevant properties of such a compact object. A localization scheme to generate 
localized excitons in clusters is proposed. Then, the local excitonic matrix elements 
as obtained in the finite clusters can be transferred to the periodic system, and, as 
the ultimate result of this approach, the excitonic bands can be calculated. 

At the end we would like to discuss some perspectives of our approach. The 
present study opens the possibility to apply highly accurate standard quantum 
chemical method, designed to account for the correlation effects in finite molecules, 
to a wide range of crystalline insulators and semiconductors being periodic systems. 
The approach provides information on the correlated wavefunction of the system 
and gives quantitatively correct results for band energies with a controlled precision 
which can systematically be improved by increasing the computational cost. With 
such a combination of highly desired features our approach has no analogue at the 
present moment. 

In the present work the method was applied to a semiconducting polymer, i.e. to 
a linear system. However, this study was going hand-in-hand with the application of 
the same methodology to bulk diamond carried out in our laboratory. By now, that 
investigation of bulk diamond has also been finished successively demonstrating that 
our approach can equally well be applied to true three-dimensional periodic infinite 
non-conducting systems. 

As discussed in section 3.4, the main computational efforts of our approach 
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goes into the evaluation of a large number of individual increments (roughly two 
hundred) which corresponds to the calculation of several hundreds of the A^-, {N — 
1)- and (A^ + l)-electron states (using a limited number of configurations only) 
providing us with the information on the energies and the wavefunctions. Such 
an amount of information requires an automatization of the calculations and data 
processing. An appropriate general-purpose computer tool, which automatically 
generates the input files for the correlation calculations with the MOLPRO program 
package, picks up all the relevant data from the output files and also computes the 
matrix elements of the effective Hamiltonians and the size-extensivity corrections 
to the correlation energies, evaluates the individual increments, checks for missing 
increments and properly sums up the increments of the correlation correction of some 
particular quantity, have been already developed in our laboratory. Also, some code 
has been designed to assist the user in generating clusters from a given crystal. 
Nonetheless, to become a routine tool used by many researchers studying solids 
a complete automatization of our method is required. By this we mean a further 
automatization of the projection of the Wannier orbitals generated by the CRYSTAL 
code onto the clusters, a dynamic and automatic control of the convergence of the 
increments with the distance between the active bonds and (by this) a proper choice 
of suitably-sized clusters for the correlation calculations, a grouping of molecular 
bonds, the summing up of the increments to both the SCF and the correlation- 
corrected local matrix elements and finally the compiling of the bands by the use 
of these LMEs. Many programm tools to tackle these tasks for some particular 
classes of systems have already been established but most of them are not yet fully 
general-purpose codes. 

On the other hand, there is an increasing activity in the implementation of 
quantum-chemical correlation methods for extended systems (MP2, say) into ex- 
isting packages for the SCF calculations which account for the electron correlation 
directly in the periodic systems using local representations of the single-particle 
wavefunctions and exploiting the predominantly local character of the correlation 
effects. The main difference of these approaches to our one is the requirement of 
all the virtual orbitals in the periodic system being localized. Then, by some cutoff 
criteria on the mutual distance of the localized orbitals the problem is effectively 
reduced to the case of correlation calculations in clusters (or some other suitable 
finite subspace). Thus, our experience will be very useful for the developers of such 
new codes. 

Let us finally turn to our approach to excitons in crystals which was estab- 
lished in the present work. Here, a completely new direction in the use of standard 
quantum-chemical program packages is opened. Further development and an ex- 
plicit implementation of this approach seems to be a very attractive task since this 
way the optical properties of non-conducting crystals can be studied by the already 
existing highly accurate quantum-chemical methods providing quantitatively cor- 
rect results. There are many open question in this new field of research and, as we 
believe, they can be answered soon. 
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